Rotated free-slip: weight the nodal normal by the facet measure the assembly integrates over (#560)

.claude/skills/adapt-on-top-faults/SKILL.md

@@ -190,6 +190,10 @@ isotropic viscosity). Rotated strong free-slip imposes `u·n̂=0` as an ESSENTIA
 constraint in a per-node (n,t) frame → machine-zero leakage AND composes with TI.
 
 ```python
+# normal=None (the default) is measure-weighted and consistent with the assembly —
+# prefer it. An analytic nhat is exact for the TRUE circle but keeps a consistency
+# error against the faceted integral (#560); use it only when the constraint must
+# follow the geometry rather than the mesh.
 nhat = mesh.CoordinateSystem.unit_e_0                    # exact radial normal (annulus/sphere)
 stokes.add_rotated_freeslip_bc(0, "Upper", normal=nhat)
 stokes.add_rotated_freeslip_bc(0, "Lower", normal=nhat)

CLAUDE.md

@@ -330,8 +330,12 @@ The PETSc-based solvers are carefully optimized and validated. **NO CHANGES with
 to impose `v·n̂ = 0`:
 
 - Enforces zero wall-normal flow to **machine precision** (Nitsche / penalty leak ~1e-3).
-- Correct on **curved / tilted / deformed** boundaries — the normal is taken per node
-  (pass an analytic `normal`, e.g. `X/|X|`, for an exact normal on curved faces).
+- Correct on **curved / tilted / deformed** boundaries — the normal is taken per node,
+  measure-weighted so it matches the straight-facet integral the assembler evaluates
+  (#560). Leave `normal=None` unless the constraint must follow the TRUE surface rather
+  than the mesh; an analytic `normal` (e.g. `X/|X|`) is exact for the geometry but keeps
+  a consistency error against the faceted assembly. See
+  `docs/developer/subsystems/rotated-freeslip.md` ("Which normal to use").
 - Works **inside the nonlinear SNES** and with **geometric FMG**. It honours
   `solver.consistent_jacobian`: use `True` (consistent Newton) for smooth nonlinear
   rheologies; `"continuation"` (staged Picard→Newton) for robustness far from the

docs/developer/subsystems/rotated-freeslip.md

@@ -16,14 +16,98 @@ stokes.add_rotated_freeslip_bc(0, "Upper", normal=nhat)          # free-slip
 stokes.add_rotated_freeslip_bc(h_dot.sym[0], "Upper", normal=nhat)  # u·n̂ = field
 ```
 
-`normal=None` uses the geometric facet normal; a sympy `1×dim` matrix in
-`mesh.X` supplies an analytic normal (exact `X/|X|` on curved boundaries — the
-preferred choice there); a constant array is also accepted. The datum must be a
+`normal=None` uses the geometric facet normal (the default, and the one
+consistent with what the assembler integrates — see below); a sympy `1×dim`
+matrix in `mesh.X` supplies an analytic normal, exact for the TRUE surface
+(`X/|X|` on a spherical cap, a constant on a planar face); a constant array is
+also accepted. The datum must be a
 *scalar* (a number, an expression of `mesh.X`, or a scalar field read); on an
 enclosed boundary it must be discretely flux-free for incompressibility. A
 corner or 3D-edge node shared between rotated boundaries has no single normal
 and stays at the free-slip pinning (the datum is ignored there).
 
+### The nodal normal is measure-weighted, not a bisector
+
+A node that sits on more than one facet — a vertex in 2D, a vertex or an
+edge-midpoint in 3D — gets one nodal normal, while the assembler integrates the
+boundary term facet by facet. The two only agree when the node's normal is
+parallel to the **measure-weighted** sum `Σ_f |f| n̂_f` (edge length in 2D, face
+area in 3D), which is what the geometric path accumulates. Plain bisector
+averaging `Σ_f n̂_f` — what UW3 did before issue #560 — is right only where the
+facets are equal; on a **kinked** wall with unequal facets it leaves a residual
+`sin(Δ/2)·(|f₁|−|f₂|)/6` in the node's free tangential row (Δ = kink angle), the
+exact constant-pressure vector stops being a null vector of the constrained
+operator, and the pressure gauge goes unpinned. **Axis-aligned** walls are
+unchanged to the last bit — their facet normals have exactly 0/±1 components, so
+`Σ_f |f| n̂_f` normalises to the same floats as `Σ_f n̂_f` whatever the weights.
+A flat but *tilted* wall is not covered by that argument and can move by one ulp.
+
+The sum runs over ALL facets meeting the node, so it must be completed **across
+ranks**. Each boundary facet is labelled on exactly one rank, so a node on a
+partition seam sees only some of its facets locally; the contributions are
+summed through the DM's local↔global scatter before normalising, which is what
+makes the normal partition-independent. Two things had to go with it: the
+outward test now points away from the facet's own support cell (the mean of the
+rank's coordinates is rank-local, and would let two facets of one node cancel),
+and the node list comes from the local mesh's exterior facets rather than the
+labelled subset, because a rank can own a node whose labelled facets are all on
+neighbours.
+
+### Which way is "outward" — and the sign of σ_nn on an inner boundary
+
+The geometric normal points away from the facet's own support cell, which is the
+**domain's** outward normal on any boundary. On a concave boundary — an annulus
+or spherical-shell **inner** arc, the CMB — that points *toward* the centre of
+curvature.
+
+This changed at #560. The old rule pointed away from the mean of the mesh
+coordinates, which on an inner arc is *into* the domain. So on a concave
+boundary, through the geometric normal:
+
+| quantity | before #560 | after |
+|---|---:|---:|
+| nodal radial component, annulus `Lower` | +1.000000 | **−1.000000** |
+| `boundary_normal_traction("Lower")` | −5.233110e-02 | **+5.233110e-02** |
+
+Magnitudes are identical to every digit; only the sign moves. `dynamic_topography_field`
+is `h = −σ_nn/(Δρ g)` on top of that number, so it reverses there too, as does the
+sign of a non-zero prescribed wall-normal datum (`u·n̂ = ũ_n`: positive now means
+outflow *from the domain* on an inner arc, where before it meant inflow). Convex
+boundaries — every box wall, an outer arc, a spherical cap — are unaffected: the two
+rules agree there.
+
+An **analytic** `normal=` is applied exactly as supplied and is *not* reoriented.
+So `X/|X|` on an inner arc is inward-of-domain and gives σ_nn of the opposite sign
+to the default. That is deliberate — the override means "use exactly this
+direction", and silently flipping it would change the meaning of a user's datum —
+but it means **you must pass `-X/|X|` on an inner boundary if you want the
+domain-outward convention.** `test_1018_rotated_nodal_normal.py` pins both halves.
+
+### Which normal to use
+
+They answer different questions, and the trade is measurable. `|A z|/|A|_F` on
+an annulus (`cellSize=0.15`), `z` = the attached constant pressure:
+
+| boundary | np | geometric (default) | analytic `X/\|X\|` |
+|---|---:|---:|---:|
+| uniform arcs | 1 | **6.6e-20** | 1.1e-14 |
+| uniform arcs | 4 | **6.7e-20** | 1.1e-14 |
+| skewed (non-uniform facets) | 1 | **6.6e-20** | 3.1e-07 |
+| skewed (non-uniform facets) | 4 | **6.7e-20** | 3.1e-07 |
+
+The geometric normal is *consistent with the assembly*: it is the direction the
+straight-facet boundary integral actually sees, so the constant pressure stays a
+null vector to machine precision at every rank count. The analytic normal is
+*consistent with the geometry*: it is tangent to the true surface, which the
+faceted mesh only approximates — so the assembler and the constraint disagree by
+an amount that grows with facet non-uniformity, and #560 does not remove it (the
+analytic column is unchanged by this fix, and identical at every rank count).
+
+Prefer the default. Reach for `normal=` when the constraint must follow the true
+surface rather than the mesh — a coarse spherical shell where faceting, not the
+gauge, is the dominant error — and be aware that the pressure gauge is then only
+as good as the numbers above.
+
 Why strong rather than Nitsche/penalty: the constraint holds to machine
 precision (a penalty leaks ~1e-3, and the leak grows exactly where anisotropy
 makes the boundary condition matter), it is correct on curved/tilted/deformed

scripts/test.sh

@@ -117,9 +117,12 @@ if [ $PARALLEL_RANKS -gt 0 ]; then
     echo "Testing global statistics and parallel operations..."
     mpirun -n $PARALLEL_RANKS python -m pytest --with-mpi tests/parallel/test_075*py || status=1
 
-    # Add other parallel test categories as they're created:
-    # echo "Testing parallel solvers..."
-    # mpirun -n $PARALLEL_RANKS python -m pytest --with-mpi tests/parallel/test_10*py || status=1
+    # Parallel SOLVER tests. This line was commented out, so test_1017 and
+    # test_1062..test_1068 — the whole rotated / constrained / MG parallel set,
+    # including the partition-independence guard for the rotated nodal normal
+    # (#560) — executed at NO rank count in CI.
+    echo "Testing parallel solvers..."
+    mpirun -n $PARALLEL_RANKS python -m pytest --with-mpi tests/parallel/test_10*py || status=1
 
     # echo "Testing parallel I/O..."
     # mpirun -n $PARALLEL_RANKS python -m pytest --with-mpi tests/parallel/test_io*py || status=1

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -6138,11 +6138,16 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         boundary : str
             Boundary label to constrain.
         normal : None or sympy 1×dim Matrix or array, optional
-            Per-node outward normal source. ``None`` uses the geometric facet
-            normal (PETSc ``computeCellGeometryFVM``; works in 2D and 3D). A
-            sympy ``1×dim`` matrix supplies an analytic normal (exact
-            ``X/|X|`` on a spherical cap, a constant on a planar face) — preferred
-            on curved boundaries. A constant array is also accepted.
+            Per-node outward normal source. ``None`` (the default, and normally
+            the right choice) uses the geometric facet normal, measure-weighted
+            so that it is consistent with the straight-facet boundary integral
+            the assembler evaluates. A sympy ``1×dim`` matrix supplies an
+            analytic normal (``X/|X|`` on a spherical cap, a constant on a
+            planar face): exact for the TRUE surface, but the assembler still
+            integrates over the facets, so it keeps a consistency error that
+            grows with facet non-uniformity. Use it when the constraint must
+            follow the geometry rather than the mesh. A constant array is also
+            accepted. See ``docs/developer/subsystems/rotated-freeslip.md``.
 
         Notes
         -----

src/underworld3/systems/free_surface.py

@@ -279,6 +279,14 @@ def _normal_direction(self, coords):
         """Per-node unit vectors along the topography direction that the surface
         increment is deformed along — vertical (Cartesian) or radial (annulus /
         spherical shell); dimension-general."""
+        # TODO(BUG): this is unconditionally +radial, but since #560 the geometric
+        # constraint normal (and so the recovered σ_nn and h_∞) is the DOMAIN's
+        # outward normal, which on a CONCAVE surface — an inner arc / CMB free
+        # surface — is −radial. The relaxation would then drive the surface away
+        # from equilibrium instead of toward it. Unreachable today: every curved
+        # free surface in the repo passes an explicit normal=rhat on an OUTER
+        # boundary, and every Cartesian one is a flat Top. Needs the sign taken from
+        # the same source as h_∞ before an inner free surface is supported.
         if self._radial:
             r = np.linalg.norm(coords, axis=1)
             r[r == 0.0] = 1.0

src/underworld3/utilities/boundary_flux.py

@@ -213,6 +213,13 @@ def _node_normals(solver, boundary, normal, nodes, dm, dim, cvec, csec, v0, v1):
     """Per-node outward unit normal (only needed to project a vector reaction).
     ``normal`` is None (geometric facet normal), a sympy 1×dim Matrix (analytic,
     lambdified), or a constant (dim,) vector."""
+    # TODO(BUG): the geometric branch below is a stale copy of the pre-#560 rule.
+    # It orients against the mean of the mesh coordinates, which is rank-local (it
+    # averages only this rank's points) and points INTO the domain on a concave
+    # boundary — rotated_bc._boundary_velocity_nodes now orients away from the
+    # facet's own support cell and sums across ranks. Currently unreachable: the
+    # only caller guards it with `if normal is not None`, so the geometric branch
+    # never runs. It will be wrong the day someone wires a geometric normal in here.
     interior_ref = cvec.mean(axis=0)
     sym_fn = const = None
     if normal is not None:

src/underworld3/utilities/fault_contact.py

@@ -548,9 +548,19 @@ def _fault_pair_nodes(solver, boundary):
               dm.getLabel(plus_name).getStratumIS(value).getIndices()
               if fS <= int(p) < fE]
 
+    # Facet normals accumulated to the pair nodes, weighted by the facet MEASURE for
+    # the same reason as the wall normals in rotated_bc (#560): the assembler
+    # integrates facet by facet, so a node on two facets is only consistent when its
+    # normal is parallel to Σ_f |f| n̂_f. The fault is not protected by being an
+    # interior surface — a constant pressure cancels exactly in the MEAN rows (the two
+    # sides carry opposite outward normals) but DOUBLES in the jump rows, and the
+    # jump-tangential (slip) row is free, so the bisector leaves √2·p·sin(δ/2)·(|f₁|−|f₂|)/6
+    # there: a pressure-driven spurious slip at every kink node, plus the same lost
+    # pressure gauge. Rank-local by construction — a seam-touching fault is
+    # redistributed onto one rank before the split, so no cross-rank sum is needed.
     nacc = {}
     for f in facets:
-        _, cent, nrm = dm.computeCellGeometryFVM(f)
+        vol, cent, nrm = dm.computeCellGeometryFVM(f)
         ne = np.asarray(nrm, dtype=float)
         ne = ne / (np.linalg.norm(ne) + 1e-30)
         support = dm.getSupport(f)
@@ -559,7 +569,7 @@ def _fault_pair_nodes(solver, boundary):
             ne = -ne
         for q in (int(c) for c in dm.getTransitiveClosure(f)[0]):
             if lsec.getFieldDof(q, _VELOCITY_FIELD) > 0:
-                nacc[q] = nacc.get(q, np.zeros(dim)) + ne
+                nacc[q] = nacc.get(q, np.zeros(dim)) + float(vol) * ne
 
     override = _compiled_normal_override(solver, boundary)
     if override is not None: