Complete the boundary-normal sum across ranks, and make the partition-independence tests self-referential (#564)

docs/developer/subsystems/rotated-freeslip.md

@@ -119,6 +119,45 @@ that must morph in time (Dirichlet→Neumann ramps).
 **Implementation**: `src/underworld3/utilities/rotated_bc.py`; registration and
 dispatch in `petsc_generic_snes_solvers.pyx` (`add_rotated_freeslip_bc`).
 
+### The same rule holds for the other two free-slip paths
+
+`add_constraint_bc` (the Lagrange-multiplier free-slip) and `add_nitsche_bc` do
+not use `_boundary_velocity_nodes`. Their default constraint direction is
+`mesh.boundary_normal(boundary)`, a P1 field assembled by
+`Mesh._assemble_boundary_normal` — a second copy of the same accumulation, and
+it was left rank-local when #560/#561 fixed the rotated one. That is **#564**:
+on `Annulus(cellSize=0.12)` the worst nodal normal on the Upper arc was 3.0e-10
+from the exact radial one in serial and **5.8e-02 (3.3°) at np=2, 3 and 4** —
+one facet's normal instead of the average of two, so its size is set by the
+facet's angular span and does not shrink with more ranks. End to end that moved
+a constrained free-slip velocity by **3.4 %** between np=1 and np=2 (#495).
+
+It is fixed the same way — the weighted contributions are summed through the
+variable's own sub-DM local↔global scatter before normalising — and all three
+accumulators now take their orientation and measure from one shared
+`utilities/facet_normals.facet_measure_and_normal`, so they cannot drift apart
+again. `boundary_flux._node_normals` is the third; its geometric branch is
+unreachable today (its caller guards it with `if normal is not None`) and it
+carries a `TODO(parallel)` rather than its own reduction.
+
+Guard: `tests/parallel/test_1069_boundary_normal_parallel.py` (an analytic
+oracle on the annulus, a global-facet-sum oracle in 2-D and 3-D, corner
+preservation, a negative control, and the Nitsche end-to-end).
+
+Two caveats that are *not* the normal, recorded so they are not re-derived:
+
+* an **internal** boundary's facets have two support cells, so
+  `Mesh._assemble_boundary_normal` skips them and
+  `mesh.boundary_normal("Internal")` comes back **zero**. `rotated_bc` keeps
+  PETSc's own face normal there instead. Neither is a supported configuration;
+  do not read either as an endorsement.
+* `add_nitsche_bc`'s default `local_h=True` scales the penalty by
+  `mesh.cell_size()`, which is **partition-dependent in its own right** (it
+  comes from a kd-tree query against this rank's centroids — see the
+  `TODO(BUG)` on `Mesh._assemble_cell_size`). With `local_h=False` the Nitsche
+  annulus agrees to 3.6e-10 at np=1…4; with it, to 6.6e-03. That is a separate
+  defect from #564 and is not fixed by it.
+
 ## One solve path
 
 There is a single driver, `solve_rotated_freeslip`: a manual outer

scripts/test.sh

@@ -79,6 +79,14 @@ if [ $PARALLEL_ONLY -eq 0 ]; then
   # no batch glob and never ran in CI.
   $PYTEST tests/test_101*py tests/test_102*py || status=1
   $PYTEST tests/test_105*py || status=1
+
+  # The boundary-normal guard lives under tests/parallel/ but carries NO
+  # mpi(min_size=2) mark, because the defect it guards is present in SERIAL in
+  # 3-D as well (#564: the facet-to-DOF routing, up to 5.9 degrees on a uniform
+  # spherical shell at np=1). Run it here so the serial job covers that path —
+  # every other serial test of the default boundary normal is on a box, where
+  # flat walls make the question vacuous. It also runs in the --p N batch below.
+  $PYTEST tests/parallel/test_1069_boundary_normal_parallel.py || status=1
   # NOT yet batched (issue #504 audit): test_106*py and test_107*py contain
   # level_2/level_3 + slow + tier_b/tier_c suites (e.g. test_1064) and need
   # a triage/deselect decision before being wired into CI.
@@ -118,9 +126,10 @@ if [ $PARALLEL_RANKS -gt 0 ]; then
     mpirun -n $PARALLEL_RANKS python -m pytest --with-mpi tests/parallel/test_075*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.
+    # test_1062..test_1069 — the whole rotated / constrained / MG parallel set,
+    # including the partition-independence guards for the rotated nodal normal
+    # (#560) and the mesh boundary normal (#564) — 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
 

src/underworld3/utilities/boundary_flux.py

@@ -29,6 +29,11 @@
 import numpy as np
 from mpi4py import MPI
 
+# One rule for a boundary facet's outward normal and measure, shared with the two
+# sibling accumulators (rotated_bc._boundary_velocity_nodes,
+# Mesh._assemble_boundary_normal).
+from underworld3.utilities.facet_normals import facet_measure_and_normal
+
 
 # M_e = (area / 12) * _P1_TRIANGLE_MASS.
 _P1_TRIANGLE_MASS = np.array(
@@ -213,14 +218,15 @@ 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)
+    # The geometric branch below now takes its orientation and its measure weight from
+    # the shared `facet_measure_and_normal` (the #560/#561 rule), so it is no longer a
+    # stale copy of the pre-#560 bisector.
+    #
+    # TODO(parallel): it still accumulates over THIS RANK's labelled facets only, so a
+    # node on a partition seam would get a partial stencil — the #564 defect. It is
+    # unreachable today (the only caller guards it with `if normal is not None`), which
+    # is why it is not carrying its own cross-rank reduction: wire one in from
+    # `rotated_bc._sum_facet_normals_across_ranks` before making this branch live.
     sym_fn = const = None
     if normal is not None:
         try:
@@ -244,13 +250,10 @@ def _node_normals(solver, boundary, normal, nodes, dm, dim, cvec, csec, v0, v1):
         for f in facets:
             if not (fS <= f < fE):
                 continue
-            _, cent, nrm = dm.computeCellGeometryFVM(f)
-            ne = np.asarray(nrm, float); ne = ne / (np.linalg.norm(ne) + 1e-30)
-            if np.dot(ne, np.asarray(cent) - interior_ref) < 0:
-                ne = -ne
+            measure, ne, _exterior = facet_measure_and_normal(dm, f)
             for q in (int(c) for c in dm.getTransitiveClosure(f)[0]):
                 if q in pts:
-                    acc[q] = acc.get(q, np.zeros(dim)) + ne
+                    acc[q] = acc.get(q, np.zeros(dim)) + measure * ne
         for q, s, _c in nodes:
             nn = acc.get(q, np.zeros(dim))
             nmap[(q, s)] = nn / (np.linalg.norm(nn) + 1e-30)

src/underworld3/utilities/facet_normals.py

@@ -0,0 +1,85 @@
+"""ONE rule for the outward normal of a boundary facet.
+
+Three places in the tree build a nodal boundary normal by walking a boundary's facets,
+taking each facet's normal from ``dm.computeCellGeometryFVM`` and accumulating it onto
+the facet's closure points:
+
+  * :func:`underworld3.utilities.rotated_bc._boundary_velocity_nodes` — the rotated
+    free-slip constraint direction,
+  * :meth:`underworld3.discretisation.Mesh._assemble_boundary_normal` — the P1 normal
+    field behind ``mesh.boundary_normal()``, which ``add_constraint_bc`` and
+    ``add_nitsche_bc`` use by default,
+  * :func:`underworld3.utilities.boundary_flux._node_normals` — the projection
+    direction for a vector reaction.
+
+They were three copies of the same six lines and they drifted: #560/#561 fixed the
+orientation rule and the measure weight in the first, leaving the other two on the
+pre-#560 rule (which orients against the mean of THIS RANK's coordinates — rank-local,
+and inward on a concave boundary). This module is that rule, once.
+
+What the rule is
+----------------
+``computeCellGeometryFVM`` returns the facet's MEASURE (edge length in 2-D, face area
+in 3-D) and a normal whose sign is PETSc's own convention, not the domain's. For an
+EXTERIOR facet the domain's outward direction is "away from the one cell the facet
+belongs to", which is local geometry — no global reference point — and is correct on a
+concave boundary (an annulus or shell inner arc), where orienting away from a
+coordinate mean points INTO the domain.
+
+Only an exterior facet has "the one cell it belongs to". An INTERNAL boundary's facets
+have two support cells and ``support[0]`` is whichever the DMPlex ordering lists first,
+so flipping against it would orient neighbouring facets of the same surface oppositely
+and they would CANCEL in a measure-weighted sum. There the raw PETSc normal is returned
+unflipped and ``exterior`` is False, and the caller decides what to do about it.
+
+The measure is returned because it is the weight a nodal normal needs to be consistent
+with the assembly that integrates the boundary term facet by facet (#560); see
+``_boundary_velocity_nodes``' docstring for the derivation.
+"""
+import numpy as np
+
+
+def facet_measure_and_normal(dm, facet, orient="exterior"):
+    """``(measure, unit normal, exterior)`` for a height-1 DMPlex point.
+
+    ``measure`` is the facet's length (2-D) / area (3-D). The normal is a unit vector
+    of length ``dm.getCoordinateDim()``. ``exterior`` is True when the facet has
+    exactly one support cell.
+
+    ``orient`` selects which of the two rules above applies:
+
+    * ``"exterior"`` (default) — flip away from the support cell ONLY when there is
+      exactly one, i.e. only where "the domain's outward direction" is defined. On an
+      internal facet the normal carries PETSc's own sign and no orientation claim is
+      made; ``exterior`` is False and the caller decides what to do about it. This is
+      what a BOUNDARY accumulator wants.
+    * ``"support0"`` — flip away from ``support[0]``'s centroid unconditionally. This
+      is what an INTERNAL SURFACE accumulator wants: on a split fault every facet of
+      the surface has the same side listed first, so "away from support[0]" is the
+      ±side split, coherent along the surface, and the reason the ``exterior`` guard
+      would be wrong there rather than merely unnecessary. Do NOT use it on a surface
+      whose facets could have their support order chosen independently — neighbouring
+      facets would then orient oppositely and CANCEL in a measure-weighted sum.
+
+    Note on dimensions: the returned vector has ``dm.getCoordinateDim()`` components
+    and the orientation dot product is taken over all of them. That is the mesh's
+    ``cdim`` for every volume mesh; on a manifold mesh (``dim < cdim``) a caller that
+    slices the result to fewer components would be taking a different dot product from
+    the one used here.
+
+    Purely rank-local: it reads geometry, takes no collective, and says nothing about
+    whether this rank sees all of the node's facets. Completing a per-node SUM across
+    ranks is the caller's job — a boundary facet is labelled on exactly one rank, so a
+    node on a partition seam sees only some of its facets locally (#564).
+    """
+    if orient not in ("exterior", "support0"):
+        raise ValueError(f"orient must be 'exterior' or 'support0', got {orient!r}")
+    measure, centroid, normal = dm.computeCellGeometryFVM(facet)
+    n = np.asarray(normal, dtype=float)
+    n = n / (np.linalg.norm(n) + 1.0e-30)
+    exterior = dm.getSupportSize(facet) == 1
+    if exterior or orient == "support0":
+        _, cell_centroid, _ = dm.computeCellGeometryFVM(int(dm.getSupport(facet)[0]))
+        if np.dot(n, np.asarray(centroid) - np.asarray(cell_centroid)) < 0.0:
+            n = -n
+    return float(measure), n, exterior

src/underworld3/utilities/fault_contact.py

@@ -48,6 +48,11 @@
 
 from underworld3 import mpi
 
+# One rule for a facet's measure and normal, shared with the three wall accumulators
+# (rotated_bc._boundary_velocity_nodes, Mesh._assemble_boundary_normal,
+# boundary_flux._node_normals). This was a fourth verbatim copy.
+from underworld3.utilities.facet_normals import facet_measure_and_normal
+
 # PETSc section field id of the velocity unknown (solver field registration
 # order: velocity first) — the same convention as rotated_bc.
 _VELOCITY_FIELD = 0
@@ -558,15 +563,18 @@ def _fault_pair_nodes(solver, boundary):
     # 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.
+    # Same rule, same one place, as the three wall accumulators — see
+    # `facet_normals.facet_measure_and_normal`. `orient="support0"` is the ONE
+    # difference and it is deliberate: a fault facet is interior by construction (two
+    # support cells), so the "exterior" rule would decline to orient it at all. Every
+    # facet of a split fault lists the same side first, so "away from support[0]" IS
+    # the ±side split and is coherent along the surface. Keeping the normalise, the
+    # +1e-30 epsilon and the measure weight shared is the point: this was a fourth
+    # verbatim copy of those six lines, and copies of them drifting apart is what
+    # produced #560 and then #564.
     nacc = {}
     for f in facets:
-        vol, cent, nrm = dm.computeCellGeometryFVM(f)
-        ne = np.asarray(nrm, dtype=float)
-        ne = ne / (np.linalg.norm(ne) + 1e-30)
-        support = dm.getSupport(f)
-        _, ccent, _ = dm.computeCellGeometryFVM(int(support[0]))
-        if np.dot(ne, np.asarray(cent) - np.asarray(ccent)) < 0:
-            ne = -ne
+        vol, ne, _exterior = facet_measure_and_normal(dm, f, orient="support0")
         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)) + float(vol) * ne

src/underworld3/utilities/rotated_bc.py

@@ -35,6 +35,11 @@
 from underworld3.utilities.boundary_flux import (
     _boundary_stratum_is, _desmear, write_boundary_scalar_field)
 
+# The outward-orientation rule for a boundary facet is shared with the two sibling
+# accumulators (Mesh._assemble_boundary_normal, boundary_flux._node_normals) so the
+# three cannot drift apart again — they did, and #564 was the bill.
+from underworld3.utilities.facet_normals import facet_measure_and_normal
+
 # The multigrid option bundles are owned in ONE place and shared with the native
 # and standard custom-P routes, so the rotated velocity block cannot be
 # configured differently from the others (#468).
@@ -222,44 +227,27 @@ def coord(q):
         if not (fS <= f < fE):
             continue
         # facet outward normal, and the measure the boundary term is integrated over
-        # (1.0 on the analytic path, whose normal does not come from the facet)
+        # (1.0 on the analytic path, whose normal does not come from the facet).
+        #
+        # The orientation rule lives in ONE place, `facet_measure_and_normal`:
+        # outward = away from the one cell this boundary facet belongs to, which is
+        # the domain's outward normal on ANY boundary, convex or not, and which needs
+        # no global reference point. Read its docstring for why the obvious
+        # alternative (away from the mean of the mesh coordinates) is both rank-local
+        # and inward on a concave boundary, and for the internal-boundary guard.
+        #
+        # CONVENTION, and it is user-visible: on a concave boundary — an annulus or
+        # shell INNER arc, the CMB — this is the opposite of what UW3 produced before
+        # #560, so σ_nn and dynamic topography read there through the GEOMETRIC normal
+        # reverse sign against earlier releases. The new sign is the one the
+        # docstrings have always claimed ("outward"). An analytic ``normal=`` is used
+        # exactly as given and is NOT reoriented, so ``X/|X|`` on an inner arc points
+        # into the domain and disagrees in sign with the default — see "Which normal
+        # to use" in ``docs/developer/subsystems/rotated-freeslip.md``.
         wgt = 1.0
         if normal is None:
-            vol, cent, nrm = dm.computeCellGeometryFVM(f)
-            ne = np.asarray(nrm, dtype=float)
-            ne = ne / (np.linalg.norm(ne) + 1e-30)
-            # Outward = away from the one cell this boundary facet belongs to, which
-            # is the domain's outward normal on ANY boundary, convex or not. The
-            # obvious alternative — away from the mean of the mesh coordinates — is
-            # BOTH rank-local (each rank averages only its own points, so two facets
-            # meeting at a seam node can be oriented oppositely and then CANCEL in
-            # the cross-rank sum) and wrong on a concave boundary, where it points
-            # INTO the domain.
-            #
-            # CONVENTION, and it is user-visible: on a concave boundary — an annulus
-            # or shell INNER arc, the CMB — this is the opposite of what UW3 produced
-            # before #560, so σ_nn and dynamic topography read there through the
-            # GEOMETRIC normal reverse sign against earlier releases. The new sign is
-            # the one the docstrings have always claimed ("outward"). An analytic
-            # ``normal=`` is used exactly as given and is NOT reoriented, so
-            # ``X/|X|`` on an inner arc points into the domain and disagrees in sign
-            # with the default — see "Which normal to use" in
-            # ``docs/developer/subsystems/rotated-freeslip.md``.
-            #
-            # Only an EXTERIOR facet has "the one cell it belongs to". An internal
-            # boundary's facets have two, and `support[0]` is whichever the DMPlex
-            # ordering happens to list first — flipping against that would orient
-            # neighbouring facets of the same surface oppositely, and they would then
-            # CANCEL in the measure-weighted sum. There the raw face normal is kept:
-            # PETSc orients it from support[0] to support[1] by its own convention,
-            # which is at least coherent along the surface. Both sibling
-            # implementations guard the same way (Mesh._assemble_boundary_normal,
-            # _local_boundary_candidates below).
-            if dm.getSupportSize(f) == 1:
-                _, ccent, _ = dm.computeCellGeometryFVM(int(dm.getSupport(f)[0]))
-                if np.dot(ne, np.asarray(cent) - np.asarray(ccent)) < 0:
-                    ne = -ne
-            wgt = float(vol)
+            vol, ne, _exterior = facet_measure_and_normal(dm, f)
+            wgt = vol
         # all velocity points on this facet (closure): verts + edges(3D) + the facet
         clo = dm.getTransitiveClosure(f)[0]
         for q in (int(c) for c in clo):