diff --git a/research/shear_refinement/README.md b/research/shear_refinement/README.md
new file mode 100644
index 0000000..67bdd8e
--- /dev/null
+++ b/research/shear_refinement/README.md
@@ -0,0 +1,70 @@
+# One measured-strain refinement example
+
+Research only: one completed base-time observation, one possible growth-only
+fluid L2 enlargement, then two more base steps. L0/L1 and the Circle-containing
+L2 rectangle stay fixed. No scheduler, coarsening, packing, CUDA changes, or
+claim of accuracy, stable long-time adaptation, or optimal mesh selection.
+
+The periodic initial shear wave is evolved for two base steps. The downloaded
+postcollision G moments supply actual L1 velocity, at dx=.5 and the parent's
+physical origin. The tagger measures |S|=sqrt(2 S:S), not vorticity. We use
+s*=|S|/0.05, cutoff .001, no wall band/history/buffer. This deliberately small
+cutoff exercises the mechanism; it is not an error estimator or calibrated
+physical criterion. A high cutoff yields no change with the same observation.
+
+The eligible corridor is L1 indices [72,76) x [20,36), immediately right of the
+fluid child. Only L1 fluid leaves are valid gradient donors and owned voters;
+covered cells, holes and ghosts cannot supply gradients. L1 has no physical
+solid cells: the Circle is entirely in the fixed L2 patch. The cropped sparse
+arrays are never treated as a full periodic domain. One-sided derivatives may
+be used at missing-neighbor edges. Unknown gradients, unavailable corridor
+cells, and tags outside the corridor are reported separately. Ignored tags do
+not certify adequate resolution elsewhere. Old coverage is always retained.
+
+`grow_rectangle` knows only a rectangle, boolean requests, alignment and area
+cap. It returns an unchanged/enlarged tuple or raises ValueError before any
+solver mutation. The caller validates the entire proposed mesh with the
+existing RemapTopology and sparse-layout guards; remap_child repeats its full
+runtime checks. This example targets exactly the original mesh once, not a
+reusable adaptation loop. A rejection propagates explicitly; no partial
+proposal is installed. All solver instances are independently closed.
+
+## Run
+
+From repository root, with dependencies already available:
+
+    PYTHONPATH=src:. python -m unittest discover -s tests/unit -p test_shear_refinement_example.py -v
+
+Only on an approved idle CUDA device:
+
+    PYTHONPATH=src:. python -m research.shear_refinement.example
+    PYTHONPATH=src:.:tests/validation/multigrid python -m unittest discover -s tests/validation/multigrid -p test_shear_refinement_example_gpu.py -v
+
+The GPU gate compares automatic selection to an explicitly identical manual
+rectangle after short continuation, requires finite leaves, exact NaN masks,
+and an unchanged source. It does not replace independent numerical replay.
+GPU execution has not been performed in the implementation cloud environment.
+
+## Tagger provenance and license
+
+`boundary_shear_tagging.py` is the unchanged accepted Grok standalone v0.2
+preparation (2026-10-16), SHA256
+`12bc70ab989953501a752a2dc23e3712b6e69bfd47e67664c1771ff2dfa4da0b`.
+The delivered `README_boundary_shear_tagging_v0.2.md` identifies it as original
+code under MIT, the same license as CelerisLab (see repository LICENSE), with
+nothing copied from upstream examples. Its source ledger/review is in the
+preparation/grok_amr_review/v0.2 handoff. Vendoring keeps this example runnable
+without external preparation paths; it is not installed as production API.
+The separate reviewed level-request policy v0.2.1 is intentionally not copied:
+this example observes only L1 and can request only the existing L2 level.
+
+## Local implementation checks (2026-10-16)
+
+Seven focused example tests passed. The combined CPU suite below passed all
+60 tests, including existing topology, pinned-body and synchronized remap gates:
+
+    PYTHONPATH=src:. python -m unittest tests.unit.test_refinement tests.unit.test_refinement_three_level_multipatch tests.unit.test_pinned_body_remap tests.unit.test_runtime_child_remap tests.unit.test_synchronized_remap tests.unit.test_shear_refinement_example -q
+
+Broad unit discovery was also attempted; 14 legacy module imports failed
+because PyCUDA is absent. That broad run is not a pass. No installation or GPU
+execution was attempted. The example uses a passive Circle, with no Jet drive.
diff --git a/research/shear_refinement/boundary_shear_tagging.py b/research/shear_refinement/boundary_shear_tagging.py
new file mode 100644
index 0000000..b5415a1
--- /dev/null
+++ b/research/shear_refinement/boundary_shear_tagging.py
@@ -0,0 +1,285 @@
+"""Boundary-proximity + physically scaled shear tagging with hysteresis (CPU prototype), v0.2.
+
+Status: standalone preparation for CelerisLab AMR (baseline dev 93d1cdd). Original code, not
+copied from any project. It only *proposes* cells for refinement: it returns masks and
+diagnostics and never builds patches, regrids, applies a max_level, or touches solver state.
+Not a block/AMR integration: callers must still legalise/align proposals and run the existing
+BlockLayout2D / FixedPatch2D / prepare_fine_boundary guards.
+
+v0.2 changes (2026-10-16 review): explicit periodic axes for distance AND dilation (minimum
+image); three separate masks (physical solid / valid stencil samples / owned cells); core tags
+only on owned-eligible fluid; hold previous core where the strain is unevaluable.
+
+Masks (all bool, shape (ny, nx)):
+  * solid  - physical solid geometry. The ONLY source of walls. Never fill holes with it.
+  * valid  - cells holding a usable, finite macroscopic velocity sample (may include ghosts).
+             Derivative stencils read only valid cells. Default: ~solid. Must not overlap solid.
+             Cells that are neither solid nor valid are "absent" (sparse holes): never walls,
+             never derivative donors, never tagged; their velocity is never read.
+  * owned  - cells this caller owns and may propose (e.g. ACTIVE fine cells, owned blocks).
+             Default: all True. eligible = owned & ~solid, and must be a subset of valid.
+  Core tags exist only on eligible cells; the buffer is clipped to owned cells.
+
+Units (match CelerisLab refinement.py: coarse dx = dt = 1, level l has dx = dt = 2**-l):
+  * Lengths (origin, wall distance, L_ref, periods) are in COARSE cells.
+  * Velocities are lattice velocities; under acoustic scaling these equal coarse units.
+  * Strain rate |S| = sqrt(2 S_ij S_ij) in 1/coarse-time-step; u = (g*y, 0) gives |S| = |g|.
+  * Shear indicator s* = |S| * L_ref / U_ref: dimensionless; the physical definition is
+    level-independent (its finite-difference estimate is exactly level-invariant only for
+    affine fields; otherwise it changes by truncation error).
+  * cell_strain = |S| * dx_level / U_ref: an analogous dx-scaled STRAIN diagnostic, returned
+    but not used to tag. It is not the quantity waLBerla/AMReX threshold: those use
+    vorticity-related indicators (rigid rotation: zero strain, nonzero vorticity). We borrow
+    only their level-scaling pattern. A dx-scaled indicator is resolution-sensitive; it is not
+    an accuracy certificate and not a substitute for a hard max_level in the policy layer.
+  * Velocity must be MACROSCOPIC (moments of one consistent population phase).
+
+Tag rule per call, on eligible cells:
+  wall:  d <= band_on or (prev and d <= band_off)               band_off >= band_on (units L_ref)
+  shear: where evaluable: s* >= shear_on or (prev and s* >= shear_off)
+         where unevaluable (no valid stencil): HOLD prev (never switches a tag off)
+  core = wall | shear;  tagged = core dilated n_buffer cells (Chebyshev, periodic-aware) & owned.
+Feed ``core`` (never ``tagged``) back as ``prev_core``. prev_core outside eligible cells is
+rejected; buffered history on an all-fluid grid cannot be detected, so the contract stays
+core-only. History is only meaningful on unchanged topology: reset or remap after a regrid.
+An untagged unevaluable cell (reported in ``unevaluable``) is NOT evidence that it may coarsen.
+"""
+from dataclasses import dataclass
+from numbers import Integral, Real
+
+import numpy as np
+
+WALL, SHEAR, HELD = 1, 2, 4  # bits in TagResult.reason (HELD: shear state held, unevaluable)
+_CHUNK = 4096
+
+
+def _ro(a):
+    a = np.ascontiguousarray(a)
+    a.setflags(write=False)
+    return a
+
+
+def _num(name, v, *, positive=False, nonneg=False):
+    if isinstance(v, (bool, np.bool_)) or not isinstance(v, Real) or not np.isfinite(v):
+        raise ValueError(f"{name} must be a finite real number")
+    v = float(v)
+    if positive and v <= 0:
+        raise ValueError(f"{name} must be > 0")
+    if nonneg and v < 0:
+        raise ValueError(f"{name} must be >= 0")
+    return v
+
+
+def _periodic(periodic):
+    if (not isinstance(periodic, (tuple, list)) or len(periodic) != 2
+            or not all(isinstance(p, (bool, np.bool_)) for p in periodic)):
+        raise ValueError("periodic must be two bools (x, y)")
+    return bool(periodic[0]), bool(periodic[1])
+
+
+def _mask(name, m, shape, default):
+    if m is None:
+        return default
+    m = np.asarray(m)
+    if m.dtype != np.bool_ or m.shape != shape:
+        raise ValueError(f"{name} must be a bool array of the grid shape")
+    return m
+
+
+def _neighbour_or(m, periodic):
+    """OR over the 3x3 Chebyshev neighbourhood; wraps on periodic axes, False-pads others."""
+    px, py = periodic
+    p = np.pad(m, ((0 if py else 1,) * 2, (0 if px else 1,) * 2), constant_values=False)
+    out = np.zeros(m.shape, bool)
+    ny, nx = m.shape
+    for dy in (-1, 0, 1):
+        for dx in (-1, 0, 1):
+            if px or py:
+                s = np.roll(p, (-dy if py else 0, -dx if px else 0), axis=(0, 1))
+            else:
+                s = p
+            oy = 0 if py else 1 + dy
+            ox = 0 if px else 1 + dx
+            out |= s[oy:oy + ny, ox:ox + nx]
+    return out
+
+
+def dilate(mask, n, periodic=(False, False)):
+    """Chebyshev dilation by n cells; wraps across periodic axes only."""
+    periodic = _periodic(periodic)
+    m = np.asarray(mask, bool).copy()
+    for _ in range(int(n)):
+        m = _neighbour_or(m, periodic)
+    return m
+
+
+@dataclass(frozen=True, eq=False)
+class TagResult:
+    core: np.ndarray          # bool: hysteresis state on eligible cells; feed back as prev_core
+    tagged: np.ndarray        # bool: core + buffer, clipped to owned; proposal for a regridder
+    reason: np.ndarray        # uint8 bits WALL|SHEAR|HELD on core cells
+    strain_star: np.ndarray   # float64 s*; NaN where unevaluable or not valid
+    cell_strain: np.ndarray   # float64 |S| dx_level / U_ref (diagnostic only)
+    wall_distance: np.ndarray # float64 coarse cells from non-solid centres; inf if no wall / solid
+    unevaluable: np.ndarray   # bool: eligible cells with no usable gradient stencil
+    eligible: np.ndarray      # bool: owned & ~solid (the only cells that can hold core tags)
+    n_unevaluable_fluid: int  # count of `unevaluable`
+
+
+def strain_rate(ux, uy, fluid, dx_level, periodic=(False, False)):
+    """|S| = sqrt(2 S:S) in 1/coarse-time-step on `fluid` (= valid sample) cells; NaN elsewhere.
+
+    Each derivative uses a central difference when both neighbours are valid (inside the
+    domain or across a periodic seam), else a one-sided difference (first order) to one valid
+    neighbour, else NaN. Values in non-valid cells are never read.
+    """
+    fluid = np.asarray(fluid, bool)
+    px, py = _periodic(periodic)
+
+    def ddir(f, axis, per):
+        n = f.shape[axis]
+        idx = np.arange(n)
+        ip, im = idx + 1, idx - 1
+        okp, okm = ip < n, im >= 0
+        if per:
+            ip, im = ip % n, im % n
+            okp = okm = np.ones(n, bool)
+        ip, im = np.clip(ip, 0, n - 1), np.clip(im, 0, n - 1)
+        take = lambda a, i: np.take(a, i, axis=axis)
+        shape = [1, 1]; shape[axis] = n
+        fp = take(fluid, ip) & okp.reshape(shape)
+        fm = take(fluid, im) & okm.reshape(shape)
+        fz = np.where(fluid, f, 0.0)
+        vp, vm = take(fz, ip), take(fz, im)
+        out = np.full(f.shape, np.nan)
+        c = fluid & fp & fm
+        out[c] = 0.5 * (vp - vm)[c]
+        o = fluid & fp & ~fm
+        out[o] = (vp - fz)[o]
+        o = fluid & ~fp & fm
+        out[o] = (fz - vm)[o]
+        return out / dx_level
+
+    dudx, dudy = ddir(ux, 1, px), ddir(ux, 0, py)
+    dvdx, dvdy = ddir(uy, 1, px), ddir(uy, 0, py)
+    sxy = 0.5 * (dudy + dvdx)
+    return np.sqrt(2.0 * (dudx**2 + dvdy**2 + 2.0 * sxy**2))
+
+
+def wall_distance(fluid, dx_level, origin=(0.0, 0.0), distance_fn=None,
+                  periodic=(False, False), solid=None):
+    """Distance (coarse cells) from each non-solid cell centre to the wall; inf if no wall.
+
+    `solid` (default ~fluid) is the physical geometry; only solid cells with a non-solid
+    Chebyshev neighbour (seam-aware) are wall cells, so absent cells are never walls.
+    distance_fn(x, y) -> signed distance (coarse units, negative inside bodies) is preferred
+    when geometry is known; it is evaluated at cell centres as given, so for periodic axes it
+    must itself return the minimum-image distance. Otherwise a raster HEURISTIC is used:
+    distance to the nearest wall-cell centre (minimum image on periodic axes) minus dx/2.
+    It is exact only for a planar wall lying on cell faces; it has no general error bound
+    (a physical wall off the cell face, curved, sub-cell or off-patch bodies differ).
+    """
+    fluid = np.asarray(fluid, bool)
+    solid = ~fluid if solid is None else np.asarray(solid, bool)
+    px, py = _periodic(periodic)
+    ny, nx = fluid.shape
+    x = origin[0] + (np.arange(nx) + 0.5) * dx_level
+    y = origin[1] + (np.arange(ny) + 0.5) * dx_level
+    X, Y = np.meshgrid(x, y)
+    target = ~solid
+    d = np.full(fluid.shape, np.inf)
+    if distance_fn is not None:
+        dv = np.asarray(distance_fn(X, Y), dtype=np.float64)
+        if dv.shape != fluid.shape or not np.all(np.isfinite(dv[target])):
+            raise ValueError("distance_fn must return finite values of the grid shape")
+        d[target] = np.maximum(dv[target], 0.0)
+        return d
+    wall = solid & _neighbour_or(~solid, (px, py))
+    if not wall.any() or not target.any():
+        return d
+    Lx, Ly = nx * dx_level, ny * dx_level
+    wx, wy = X[wall], Y[wall]
+    fx, fy = X[target], Y[target]
+    best = np.empty(fx.size)
+    for s in range(0, fx.size, _CHUNK):
+        ex, ey = np.abs(fx[s:s + _CHUNK, None] - wx), np.abs(fy[s:s + _CHUNK, None] - wy)
+        if px:
+            ex = np.minimum(ex, Lx - ex)
+        if py:
+            ey = np.minimum(ey, Ly - ey)
+        best[s:s + _CHUNK] = np.sqrt((ex * ex + ey * ey).min(axis=1))
+    d[target] = np.maximum(best - 0.5 * dx_level, 0.0)
+    return d
+
+
+def tag_cells(ux, uy, solid, *, dx_level, U_ref, L_ref, shear_on, shear_off,
+              band_on, band_off, n_buffer=2, prev_core=None, origin=(0.0, 0.0),
+              periodic=(False, False), distance_fn=None, valid=None, owned=None):
+    """Propose refinement tags for one level. Pure function: inputs are not modified.
+
+    band_on/band_off are in units of L_ref. See module docstring for the three masks.
+    """
+    solid = np.asarray(solid)
+    if solid.dtype != np.bool_ or solid.ndim != 2 or min(solid.shape) < 3:
+        raise ValueError("solid must be a 2-D bool array with both sides >= 3")
+    shape = solid.shape
+    ux = np.asarray(ux, dtype=np.float64)
+    uy = np.asarray(uy, dtype=np.float64)
+    if ux.shape != shape or uy.shape != shape:
+        raise ValueError("ux, uy and solid must have the same (ny, nx) shape")
+    valid = _mask("valid", valid, shape, ~solid)
+    owned = _mask("owned", owned, shape, np.ones(shape, bool))
+    if np.any(valid & solid):
+        raise ValueError("valid must not overlap solid (solid cells are not velocity samples)")
+    eligible = owned & ~solid
+    if np.any(eligible & ~valid):
+        raise ValueError("owned non-solid cells must be valid (have a velocity sample)")
+    if not (np.all(np.isfinite(ux[valid])) and np.all(np.isfinite(uy[valid]))):
+        raise ValueError("velocity must be finite on every valid cell")
+    dx_level = _num("dx_level", dx_level, positive=True)
+    lvl = -np.log2(dx_level)
+    if abs(lvl - round(lvl)) > 1e-12 or round(lvl) < 0:
+        raise ValueError("dx_level must be 2**-level coarse cells (1, 0.5, 0.25, ...)")
+    U_ref = _num("U_ref", U_ref, positive=True)
+    L_ref = _num("L_ref", L_ref, positive=True)
+    shear_on = _num("shear_on", shear_on, positive=True)
+    shear_off = _num("shear_off", shear_off, nonneg=True)
+    if not shear_off < shear_on:
+        raise ValueError("hysteresis needs shear_off < shear_on")
+    band_on = _num("band_on", band_on, nonneg=True)
+    band_off = _num("band_off", band_off, nonneg=True)
+    if not band_off >= band_on:
+        raise ValueError("hysteresis needs band_off >= band_on")
+    if isinstance(n_buffer, (bool, np.bool_)) or not isinstance(n_buffer, Integral) or n_buffer < 0:
+        raise ValueError("n_buffer must be a non-negative integer")
+    if len(origin) != 2:
+        raise ValueError("origin must be (x0, y0) in coarse units")
+    origin = (_num("origin[0]", origin[0]), _num("origin[1]", origin[1]))
+    periodic = _periodic(periodic)
+    if prev_core is None:
+        prev = np.zeros(shape, bool)
+    else:
+        prev = np.asarray(prev_core)
+        if prev.dtype != np.bool_ or prev.shape != shape:
+            raise ValueError("prev_core must be a bool array of the grid shape")
+        if np.any(prev & ~eligible):
+            raise ValueError("prev_core marks solid/unowned cells: pass TagResult.core, not tagged")
+
+    S = strain_rate(ux, uy, valid, dx_level, periodic)
+    s_star = S * L_ref / U_ref
+    cell_strain = S * dx_level / U_ref
+    d = wall_distance(~solid, dx_level, origin, distance_fn, periodic, solid=solid)
+    dn = d / L_ref
+
+    ok = np.isfinite(s_star)
+    sv = np.where(ok, s_star, -np.inf)
+    unevaluable = eligible & ~ok
+    held = unevaluable & prev
+    shear = eligible & (((sv >= shear_on) | (prev & (sv >= shear_off))) | held)
+    wall = eligible & ((dn <= band_on) | (prev & (dn <= band_off)))
+    core = shear | wall
+    reason = (wall * WALL + shear * SHEAR + held * HELD).astype(np.uint8)
+    tagged = dilate(core, n_buffer, periodic) & owned
+    return TagResult(_ro(core), _ro(tagged), _ro(reason), _ro(s_star), _ro(cell_strain),
+                     _ro(d), _ro(unevaluable), _ro(eligible),
+                     int(np.count_nonzero(unevaluable)))
diff --git a/research/shear_refinement/example.py b/research/shear_refinement/example.py
new file mode 100644
index 0000000..66295d3
--- /dev/null
+++ b/research/shear_refinement/example.py
@@ -0,0 +1,140 @@
+"""One measured L1 strain observation, one growth-only L2 remap, no scheduler.
+
+Run from the repository root: PYTHONPATH=src:. python -m research.shear_refinement.example
+Requires CUDA; CPU helpers and tests never import CUDA.
+"""
+from contextlib import ExitStack
+from numbers import Integral
+import numpy as np
+
+from CelerisLab.config import LBMConfig
+from CelerisLab.body.geometry.circle import CircleGeometry
+from CelerisLab.body.objects import SimObject
+from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
+from CelerisLab.lbm.refinement import FixedPatchSet2D
+from CelerisLab.lbm.synchronized_remap import RemapTopology
+from .boundary_shear_tagging import tag_cells
+
+CFG = LBMConfig(nx=48, ny=32, viscosity=.1)
+PARENTS = ((4, 4, 20, 28), (24, 4, 44, 28))
+CHILDREN = ((16, 16, 32, 40), (56, 16, 72, 40))
+CORRIDOR = (72, 20, 76, 36)  # L1 allocated indices; bounded to right of fluid L2
+
+
+def grow_rectangle(current, request, *, alignment=4, area_cap=480):
+    """Bounding-box growth only, in request-array cell-edge coordinates.
+
+    No body or solver knowledge. Reject invalid/budget-exceeding proposals with
+    ValueError; never clip a request, shrink old coverage or mutate an input.
+    Alignment is edge divisibility, not a substitute for full topology guards.
+    """
+    request = np.asarray(request)
+    if request.dtype != bool or request.ndim != 2:
+        raise ValueError('request must be a 2D boolean mask')
+    values = (*current, alignment, area_cap)
+    if len(current) != 4 or any(isinstance(v, (bool, np.bool_)) or not isinstance(v, Integral) for v in values):
+        raise ValueError('rectangle, alignment and area cap must be integers')
+    x0, y0, x1, y1 = current
+    ny, nx = request.shape
+    if alignment <= 0 or area_cap <= 0 or not (0 <= x0 < x1 <= nx and 0 <= y0 < y1 <= ny):
+        raise ValueError('invalid rectangle, alignment or area cap')
+    if any(v % alignment for v in current):
+        raise ValueError('current rectangle must be aligned')
+    y, x = np.nonzero(request)
+    if x.size:
+        x0, y0 = min(x0, int(x.min())), min(y0, int(y.min()))
+        x1, y1 = max(x1, int(x.max())+1), max(y1, int(y.max())+1)
+    a = alignment
+    result = (x0//a*a, y0//a*a, (x1+a-1)//a*a, (y1+a-1)//a*a)
+    if result[2] > nx or result[3] > ny:
+        raise ValueError('aligned growth exceeds observation frame')
+    if (result[2]-result[0])*(result[3]-result[1]) > area_cap:
+        raise ValueError('growth exceeds area cap')
+    return result
+
+
+def objects():
+    return (SimObject(71, geometry=CircleGeometry(12., 14., .8),
+                      center=(12., 14.), radius=.8),)
+
+
+def topology(rectangles=CHILDREN):
+    return RemapTopology.from_config(CFG, FixedPatchSet2D.from_config(CFG, PARENTS),
+                                     rectangles, child_objects=objects())
+
+
+def observe(g, mesh, *, shear_on=.001):
+    """Moments of completed-time postcollision G; only real L1 fluid leaves.
+
+    No reconstructed covered values, ghosts or holes are derivative donors.
+    solid is independent geometry (empty on L1: the body is pinned at L2).
+    The array is cropped/sparse, so neither derivative axis wraps.
+    """
+    valid = mesh.leaf_masks[1]
+    f = g[:, valid].astype(np.float64)
+    rho = f.sum(axis=0)
+    if not np.isfinite(f).all() or np.any(rho <= 0):
+        raise ValueError('invalid completed-time L1 fluid populations')
+    ux, uy = (np.full(valid.shape, np.nan) for _ in range(2))
+    ux[valid] = (np.asarray(D2Q9_EX)[:, None]*f).sum(axis=0)/rho
+    uy[valid] = (np.asarray(D2Q9_EY)[:, None]*f).sum(axis=0)/rho
+    tags = tag_cells(ux, uy, np.zeros(valid.shape, bool), valid=valid, owned=valid,
+                     dx_level=.5, origin=mesh.parent.fine_origin, periodic=(False, False),
+                     U_ref=.05, L_ref=1., shear_on=shear_on, shear_off=0.,
+                     band_on=0., band_off=0., n_buffer=0)
+    corridor = np.zeros(valid.shape, bool)
+    x0, y0, x1, y1 = CORRIDOR
+    corridor[y0:y1, x0:x1] = True
+    request = tags.tagged & corridor
+    proposed = (CHILDREN[0], grow_rectangle(CHILDREN[1], request))
+    # Host guards cover nesting, donor footprints, complete body and separation.
+    candidate = topology(proposed)
+    candidate.child.fine_layout((4, 4))
+    report = dict(requested=int(request.sum()), unknown=int(tags.unevaluable.sum()),
+                  corridor_unknown=int((tags.unevaluable & corridor).sum()),
+                  corridor_unavailable=int((corridor & ~valid).sum()),
+                  out_of_scope=int((tags.tagged & ~corridor).sum()), rectangles=proposed)
+    return proposed, report
+
+
+def initial_fields(mesh):
+    """Smooth periodic shear wave, sampled in each level's physical frame."""
+    from CelerisLab.lbm.equilibrium import compute_feq_d2q9
+    frames = ((np.zeros(2), 1.), (np.asarray(mesh.parent.fine_origin), .5),
+              (np.asarray(mesh.parent.fine_origin)+.5*np.asarray(mesh.child.fine_origin), .25))
+    arrays = []
+    for leaf, (origin, dx) in zip(mesh.leaf_masks, frames):
+        y = origin[1]+dx*(np.arange(leaf.shape[0])+.5)
+        ux = .01*np.sin(2*np.pi*y/CFG.ny)
+        rows = np.stack([compute_feq_d2q9(1., u, 0.) for u in ux], axis=1)
+        arrays.append(np.broadcast_to(rows[:, :, None], (9, *leaf.shape)).copy())
+    return tuple(arrays)
+
+
+def initialize(run):
+    run.initialize(*initial_fields(topology()))
+
+
+def make(context):
+    from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+    return FixedPatchStepper(CFG, PARENTS, context=context, base_block=(4, 4),
+                             fine_block=(4, 4), child_block=(4, 4),
+                             child_rectangles_in_parent=CHILDREN, child_objects=objects())
+
+
+def main():
+    from CelerisLab.cuda.context import CudaContext
+    with CudaContext(0) as ctx, ExitStack() as owned:
+        run = owned.enter_context(make(ctx))
+        initialize(run)
+        run.step(2)  # Returns only at completed synchronized base time.
+        proposed, report = observe(run.download()[1], topology())
+        print(report)
+        if proposed != CHILDREN:
+            run = owned.enter_context(run.remap_child(proposed, child_objects=objects()))
+        run.step(2)
+        print(f'completed base step {run.step_count}; body rectangle {run.child_patch.bounds[0]}')
+
+
+if __name__ == '__main__':
+    main()
diff --git a/tests/unit/test_shear_refinement_example.py b/tests/unit/test_shear_refinement_example.py
new file mode 100644
index 0000000..a96d173
--- /dev/null
+++ b/tests/unit/test_shear_refinement_example.py
@@ -0,0 +1,84 @@
+"""CPU-only research-example contracts; no solver/GPU import required."""
+import unittest
+import numpy as np
+from research.shear_refinement.boundary_shear_tagging import tag_cells, strain_rate
+from research.shear_refinement.example import (CHILDREN, grow_rectangle, topology,
+                                                observe, initial_fields)
+
+
+class ShearRefinementExample(unittest.TestCase):
+    def test_affine_shear_physical_spacing_and_rotation(self):
+        y, x = np.indices((12, 20), dtype=float)
+        x, y = 7.+.5*(x+.5), 3.+.5*(y+.5)
+        valid = np.ones(x.shape, bool)
+        np.testing.assert_allclose(strain_rate(.003*y, 0*x, valid, .5), .003, atol=1e-17)
+        np.testing.assert_allclose(strain_rate(-.003*y, .003*x, valid, .5), 0., atol=1e-17)
+
+    def test_holes_are_unknown_not_walls_or_zero(self):
+        valid = np.zeros((7, 7), bool); valid[3, 3] = True
+        u = np.full(valid.shape, np.nan); u[3, 3] = .1
+        result = tag_cells(u, u, np.zeros_like(valid), valid=valid, owned=valid,
+                           dx_level=.5, origin=(7., 3.), U_ref=.05, L_ref=1.,
+                           shear_on=.01, shear_off=0., band_on=0., band_off=0., n_buffer=0)
+        self.assertTrue(result.unevaluable[3, 3])
+        self.assertTrue(np.isnan(result.strain_star[3, 3]))
+        self.assertFalse(result.tagged.any())
+        self.assertTrue(np.isinf(result.wall_distance).all())
+
+    def test_masked_hole_never_supplies_a_gradient(self):
+        y, x = np.indices((9, 9), dtype=float)
+        valid = np.ones(x.shape, bool); valid[:, 4] = False
+        ux = .003*(5.+.5*y); uy = np.zeros_like(ux)
+        ux[~valid] = np.inf; uy[~valid] = np.nan
+        strain = strain_rate(ux, uy, valid, .5)
+        np.testing.assert_allclose(strain[valid], .003, atol=1e-17)
+        self.assertTrue(np.isnan(strain[~valid]).all())
+
+    def test_predictable_growth_and_input_unchanged(self):
+        mask = np.zeros((64, 96), bool); mask[22, 74] = True
+        before = mask.copy()
+        self.assertEqual(grow_rectangle(CHILDREN[1], mask), (56, 16, 76, 40))
+        np.testing.assert_array_equal(mask, before)
+        mask[:] = False
+        self.assertEqual(grow_rectangle(CHILDREN[1], mask), CHILDREN[1])
+
+    def test_budget_and_geometry_rejection(self):
+        mask = np.zeros((64, 96), bool); mask[22, 74] = True
+        with self.assertRaisesRegex(ValueError, 'area cap'):
+            grow_rectangle(CHILDREN[1], mask, area_cap=479)
+        with self.assertRaises(ValueError):
+            grow_rectangle((57, 16, 72, 40), mask)
+        mask = np.zeros((9, 9), bool); mask[8, 8] = True
+        with self.assertRaisesRegex(ValueError, 'frame'):
+            grow_rectangle((0, 0, 4, 4), mask)
+        with self.assertRaises(ValueError):
+            topology((CHILDREN[0], (56, 16, 88, 40)))
+
+    def test_uniform_velocity_observation_does_not_grow(self):
+        from CelerisLab.lbm.equilibrium import compute_feq_d2q9
+        mesh = topology()
+        g = np.broadcast_to(compute_feq_d2q9(1., .01, 0.)[:, None, None],
+                            (9, *mesh.leaf_masks[1].shape)).copy()
+        g[:, ~mesh.leaf_masks[1]] = np.nan
+        rectangles, report = observe(g, mesh)
+        self.assertEqual(rectangles, CHILDREN)
+        self.assertEqual(report['requested'], 0)
+
+    def test_observation_growth_and_low_strain_no_change(self):
+        mesh = topology()
+        g = initial_fields(mesh)[1]
+        g[:, ~mesh.leaf_masks[1]] = np.nan
+        before = g.copy()
+        rectangles, report = observe(g, mesh)
+        self.assertEqual(rectangles, (CHILDREN[0], (56, 16, 76, 40)))
+        self.assertEqual(report['requested'], 64)
+        self.assertEqual(report['corridor_unknown'], 0)
+        self.assertEqual(report['corridor_unavailable'], 0)
+        self.assertGreater(report['out_of_scope'], 0)
+        np.testing.assert_array_equal(g, before)
+        rectangles, report = observe(g, mesh, shear_on=1.)
+        self.assertEqual(rectangles, CHILDREN)
+        self.assertEqual(report['requested'], 0)
+
+
+if __name__ == '__main__': unittest.main()
diff --git a/tests/validation/multigrid/test_shear_refinement_example_gpu.py b/tests/validation/multigrid/test_shear_refinement_example_gpu.py
new file mode 100644
index 0000000..a02fc71
--- /dev/null
+++ b/tests/validation/multigrid/test_shear_refinement_example_gpu.py
@@ -0,0 +1,43 @@
+"""Run explicitly on an idle GPU; automatic/manual identical-layout equivalence.
+
+This is not independent flow-accuracy validation. CPU collection never runs it.
+PYTHONPATH=src:.:tests/validation/multigrid python -m unittest discover \
+    -s tests/validation/multigrid -p test_shear_refinement_example_gpu.py -v
+"""
+from contextlib import ExitStack
+import unittest
+import numpy as np
+from CelerisLab.cuda.context import CudaContext
+from research.shear_refinement.example import (CHILDREN, make, initialize, observe,
+                                                topology, objects)
+
+
+class ShearRefinementGPU(unittest.TestCase):
+    def test_measured_growth_matches_explicit_manual_continuation(self):
+        with CudaContext(0) as ctx, ExitStack() as owned:
+            old = owned.enter_context(make(ctx))
+            initialize(old)
+            old.step(2)
+            before = old.download()
+            chosen, report = observe(before[1], topology())
+            manual = (CHILDREN[0], (56, 16, 76, 40))
+            self.assertEqual(chosen, manual)
+            self.assertGreater(report['requested'], 0)
+            self.assertEqual(report['corridor_unknown'], 0)
+            self.assertEqual(report['corridor_unavailable'], 0)
+            auto = owned.enter_context(old.remap_child(chosen, child_objects=objects()))
+            explicit = owned.enter_context(old.remap_child(manual, child_objects=objects()))
+            for run in (auto, explicit):
+                self.assertEqual(run.child_patch.bounds[0], CHILDREN[0])
+                self.assertEqual(run.step_count, 2)
+                run.step(2)
+            for a, b, leaf in zip(auto.download(), explicit.download(), topology(manual).leaf_masks):
+                np.testing.assert_array_equal(a, b)
+                self.assertTrue(np.isfinite(a[:, leaf]).all())
+                self.assertTrue(np.isnan(a[:, ~leaf]).all())
+            for a, b in zip(old.download(), before):
+                np.testing.assert_array_equal(a, b)
+            self.assertEqual((old.step_count, auto.step_count, explicit.step_count), (2, 4, 4))
+
+
+if __name__ == '__main__': unittest.main()