diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md index f1ed4cc13..83606e1ca 100644 --- a/docs/developer/subsystems/rotated-freeslip.md +++ b/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 diff --git a/scripts/test.sh b/scripts/test.sh index 569ab4558..0f70ae170 100755 --- a/scripts/test.sh +++ b/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 diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 03fe527f2..ecaf3eb1c 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2909,7 +2909,7 @@ def _update_projected_normals(self): def boundary_normal(self, boundary): """Outward unit normal of a single boundary, tracking deformation. - Assembles the EXACT, outward, area-weighted PETSc facet normals + Assembles the EXACT, outward, measure-weighted PETSc facet normals (``dm.computeCellGeometryFVM``) from ONLY this boundary's facets onto its P1 vertices. Because each boundary is assembled independently, a vertex shared by two boundaries (a sharp corner) is NOT averaged @@ -2917,6 +2917,31 @@ def boundary_normal(self, boundary): smooth boundary (e.g. a free surface) the result is the smooth deformed normal. Cached per boundary; rebuilt lazily after a deform. + COLLECTIVE. The per-vertex sum runs over ALL the facets meeting the + vertex, which in parallel are split between ranks, so it is completed + through the variable's own local↔global scatter before normalising + (#564). Every rank must call this, including one that owns no part of + the boundary. + + .. note:: **The 3-D answer changed at #564, in SERIAL as well as in + parallel.** Facet contributions used to be routed to "the DOFs + nearest the facet centroid" by a kd-tree query. On a TETRAHEDRAL + boundary the three DOFs nearest a face centroid are not always that + face's own three vertices, so the query silently picked up a + neighbour and the assembled normal was wrong — by up to **0.103 in + the unit normal (≈5.9°) on ``SphericalShell(0.55, 1.0, cs=0.35)``, + at np=1, on a UNIFORM mesh**. Rows now come from the variable's own + section, which is exact: measured against an independent + global-facet-sum oracle, 1.9e-16 (Upper) and 2.2e-16 (Lower) where the + old route was 4.7e-02 and 1.03e-01. + + 2-D is unaffected (annulus 1.1e-16 old and new, box bit-identical) — + on an edge the two nearest DOFs to the midpoint are always its own two + vertices. So any **3-D** result that used the default ``normal=None`` + on a curved boundary moves, and moves toward the right answer. Guarded + by ``tests/parallel/test_1069_boundary_normal_parallel.py``, which runs + at np=1 as well as np>1 precisely so this path is covered. + Returns a sympy Matrix (row) of the P1 normal-field components, for use as the constraint direction in Nitsche/penalty BCs. @@ -2975,56 +3000,182 @@ def _boundary_facets(self, name): return face_pts def _assemble_boundary_normal(self, var, name): - """Fill ``var`` with the area-weighted outward facet normal assembled - from the faces of boundary ``name`` only (see :meth:`boundary_normal`).""" - from scipy.spatial import cKDTree + """Fill ``var`` with the measure-weighted outward facet normal assembled + from the faces of boundary ``name`` only (see :meth:`boundary_normal`). + + The nodal normal is ``Σ_f |f| n̂_f`` over EVERY facet of this boundary that + meets the node, normalised at the end. In parallel that sum has to be + completed across ranks before the normalise: a boundary facet is labelled + on exactly one rank (measured — see the note below), so a node on a + partition seam sees only SOME of its facets locally and normalising a + partial stencil gives it a rotated normal. That is #564: on an + ``Annulus(cellSize=0.12)`` the Upper arc's worst nodal normal was 3.0e-10 + from the exact radial one in serial and 5.8e-02 (3.3 degrees) at np=2,3,4 — + exactly the error of taking one facet's normal instead of the average of + the two — and it moved a constrained free-slip answer by 3.4 %. + + COLLECTIVE. The reduction runs on every rank, including one that owns no + facet of this boundary (the #405 lesson: a rank-local early return here + deadlocks the ranks that do own facets). + + Why a plain ADD is exact, with no de-duplication: no boundary facet is + labelled on two ranks and none is labelled away from its owner. Measured + on the annulus and the spherical shell at np=2,3,4 for BOTH label sources + (the per-boundary label this uses and the consolidated ``UW_Boundaries``); + the guard that keeps it true is + ``tests/parallel/test_1069_boundary_normal_parallel.py``. + + FAILURE IS COLLECTIVE AND LOUD. A rank that cannot complete its facet walk + raises :class:`RuntimeError` on EVERY rank, not just its own. Two failure + modes are being avoided, and both were measured on the first version of + this routine: + + * swallowing the failure rank-locally and carrying on gives that rank's + OWNED boundary DOFs a ZERO normal — so the constraint direction over + that part of the boundary is the zero vector, with a converged solve + and no message. That is strictly worse than the 3.3-degree error this + routine exists to remove, and indistinguishable from success; + * raising rank-locally takes that rank out of the sub-DM collectives + below and HANGS the others (measured: rank 1 returns, rank 0 blocks, + killed by the launcher timeout). + + So the flag is agreed with an all-reduce first, every rank then takes the + same branch, and every rank raises. Callers that swallow the exception + (:meth:`deform`) are therefore safe by construction: what reaches them is + already symmetric. + """ + from underworld3.utilities.facet_normals import facet_measure_and_normal + cdim = self.cdim dm = self.dm - coords = numpy.ascontiguousarray(var.coords) - accum = numpy.zeros((coords.shape[0], cdim)) + comm = dm.comm.tompi4py() + failed = 0 + detail = "" + ncomp = None + accum = None + + # Rank-local set-up. Guarded like the facet walk below and for the same + # reason: `dm.createSubDM` is COLLECTIVE, so a rank must not leave before + # reaching it. + try: + ncomp = var.num_components + # One node per DMPlex point is assumed by the `offset // ncomp` row + # arithmetic below. True for the degree-1 variable `boundary_normal` + # builds, but that path adopts a pre-existing `_n_bd_` variable + # if one is already registered (a checkpoint restore, or user code), + # and a higher-degree space puts several nodes on one point — a P3 + # edge carries two in 2-D — which would land both on the first row and + # leave the second at zero. Refuse rather than silently half-fill. + if var.degree != 1: + raise RuntimeError( + f"boundary_normal needs a degree-1 field; '{var.clean_name}' is " + f"degree {var.degree}. A higher-degree space carries several " + f"nodes per DMPlex point and this assembly writes one row per " + f"point.") + # Dense over this rank's LOCAL DOFs (ghosts included), so a node whose + # labelled facets all live on a neighbour needs no special enumeration — + # it is simply a row that stays zero until the reduction fills it. + accum = numpy.zeros_like(numpy.asarray(var.data)) + except Exception as exc: + failed, detail = 1, f"{type(exc).__name__}: {exc}" + + # DOF rows come from the variable's OWN section on its sub-DM: exact on + # every rank, and the same section the reduction below scatters through. + # (This used to be a kd-tree lookup of the DOFs nearest the facet + # centroid. That is a heuristic, and not only on a graded mesh: on a + # TETRAHEDRAL boundary the three DOFs nearest a face centroid are not + # always that face's own three vertices, so it mis-assigned on a uniform + # spherical shell in SERIAL — see :meth:`boundary_normal`. It also cannot + # be made to agree across ranks, because each rank's tree is built from + # its own local coordinates.) + indexset, subdm = dm.createSubDM(var.field_id) + try: + try: + ssec = subdm.getLocalSection() + for f in self._boundary_facets(name): + # Orientation needs the facet's OWN support cell, and only an + # exterior facet has exactly one. An internal boundary's facets + # have two and support[0] is arbitrary, so neighbouring facets + # of the same surface could be oriented oppositely and CANCEL + # in the sum — those facets are skipped, which is why a normal + # requested for an INTERNAL boundary comes back zero. + # (rotated_bc keeps the raw PETSc normal there instead; + # neither is a supported use of this routine today.) + measure, nrm, exterior = facet_measure_and_normal(dm, f) + if not exterior: + continue + for q in (int(c) for c in dm.getTransitiveClosure(f)[0]): + if ssec.getDof(q) <= 0: + continue + accum[ssec.getOffset(q) // ncomp] += measure * nrm[:cdim] + except Exception as exc: + if not failed: + failed, detail = 1, f"{type(exc).__name__}: {exc}" + + # Agree on the outcome BEFORE the reduction, so every rank takes the + # same branch. Skipping the reduction on a failure is what keeps the + # raise below from being reached by only some ranks. + if comm.size > 1: + failed = comm.allreduce(failed) + if not failed: + accum = self._sum_local_dofs_across_ranks(subdm, accum) \ + if comm.size > 1 else accum + finally: + indexset.destroy() + subdm.destroy() - face_pts = self._boundary_facets(name) + if failed: + reports = [d for d in (comm.allgather(detail) if comm.size > 1 + else [detail]) if d] + raise RuntimeError( + f"boundary normal assembly for {name!r} failed on " + f"{failed} of {comm.size} rank(s): {'; '.join(reports[:4])}") - tree = cKDTree(coords) - # P1 vertices per facet, counted from the facet's own closure so this - # works for non-simplex facets too (2D edge=2, 3D tri=3, 3D quad=4). - vStart, vEnd = dm.getDepthStratum(0) - for f in face_pts: - if dm.getSupportSize(f) != 1: - continue - area, cent, nrm = dm.computeCellGeometryFVM(f) - nrm = numpy.asarray(nrm)[:cdim] - cell = dm.getSupport(f)[0] - _, ccent, _ = dm.computeCellGeometryFVM(cell) - if numpy.dot(nrm, numpy.asarray(cent)[:cdim] - - numpy.asarray(ccent)[:cdim]) < 0: - nrm = -nrm - _clo = dm.getTransitiveClosure(f)[0] - nverts = int(numpy.count_nonzero((_clo >= vStart) & (_clo < vEnd))) or cdim - # Accumulate to the facet's P1 DOFs (its vertices) — found as the - # nearest DOFs to the facet centroid. This avoids indexing the local - # coordinate array by (vertex_point - vStart), which is only valid - # in serial (the parallel coordinate layout differs → out-of-range). - # Normalisation at the end makes the per-DOF weight (full vs share) - # irrelevant to the resulting direction. - _, idxs = tree.query(numpy.asarray(cent)[:cdim], k=nverts) - for idx in numpy.atleast_1d(idxs): - accum[idx] += area * nrm - - # TODO(parallel): a boundary vertex on a partition interface should - # ADD-reduce the UNnormalised facet contributions from both ranks before - # normalising (DMLocalToGlobal ADD_VALUES on the variable's section), - # so its normal is the full-stencil average rather than this rank's - # partial stencil. This is parallel-SAFE as-is (rank-interior boundary - # vertices are exact; only the handful of partition-seam surface - # vertices get a slightly-rotated unit normal). A first ADD-reduce - # attempt SEGV'd on the lazily-built work variable's global vec; deferred - # to a focused follow-up with the right vec/section plumbing. mag = numpy.sqrt(numpy.sum(accum ** 2, axis=1)) nonzero = mag > 1.0e-30 accum[nonzero] /= mag[nonzero, numpy.newaxis] var.data[...] = accum + def _sum_local_dofs_across_ranks(self, subdm, values): + """Complete a per-DOF sum across ranks: this rank's own contributions in, + the sum over EVERY rank's contributions out — the same value on every rank + that holds the DOF. + + ``values`` is dense over ``subdm``'s LOCAL DOFs, shaped ``(ndof, ncomp)`` + exactly like the variable's ``.data``. The sum rides the sub-DM's own + local↔global scatter: ADD into the global vector accumulates every ghost + copy onto the owner, and scattering back hands every rank the identical + total. Work vectors are created here rather than borrowed from the + variable, whose global vec is built lazily. + + COLLECTIVE on the DM's communicator. + + There is deliberately NO "the round trip lost this DOF, keep the local + value" fallback. The question such a fallback wants to ask is "was this + DOF constrained out of the global vector?"; the only thing it can cheaply + test is "did it come back all-zero?", and those two differ on exactly the + input where it would matter — a node whose GLOBAL contributions cancel + (opposed facets on a degenerate or zero-thickness boundary) would have its + rank-local PARTIAL value restored, re-introducing #564 on the one case the + reduction exists to get right. The mesh DM carries no essential-BC + constraints (those live on the solver DM), so nothing is lost today; if + that ever changes, this should fail loudly rather than guess, and + ``ssec.getConstraintDof(q)`` is the predicate to use. + """ + ncomp = values.shape[1] + lvec = subdm.createLocalVector() + gvec = subdm.createGlobalVector() + try: + lvec.array[...] = values.reshape(-1) + gvec.set(0.0) + subdm.localToGlobal(lvec, gvec, addv=PETSc.InsertMode.ADD_VALUES) + subdm.globalToLocal(gvec, lvec, addv=PETSc.InsertMode.INSERT_VALUES) + summed = numpy.array(lvec.array, dtype=float).reshape(-1, ncomp) + finally: + lvec.destroy() + gvec.destroy() + return summed + def cell_size(self): """Local, per-cell characteristic mesh size as a scalar field symbol. @@ -3106,6 +3257,23 @@ def _assemble_cell_size(self, var): access, no collective): mixing a rank-local fast path with a collective fallback would diverge across ranks and deadlock, because ``var.coords`` triggers the collective ``_get_coords_for_basis``.""" + # TODO(BUG): this field is PARTITION-DEPENDENT, and so therefore is the + # Nitsche penalty gamma*mu/h that consumes it (local_h=True, the default). + # Not the indexing here — the values. `_get_mesh_sizes` measures a cell by + # the distance from its vertices to the NEAREST CENTROID in a kd-tree built + # from THIS RANK's centroids, so near a partition seam the nearest centroid + # may simply be absent. Measured on Annulus(cellSize=0.12): the field's sum + # is 26.0822 at np=1, 26.1211 at np=2 and 26.1386 at np=4, and its max moves + # at np=4. End to end that is 6.6e-03 in the velocity of a Nitsche free-slip + # annulus and it does NOT shrink with solver tolerance. + # This is a DIFFERENT defect from the boundary normal fixed for #564 (which + # is now clean: the same solve with local_h=False agrees to 3.6e-10 at + # np=1..4). It is the local h that is left, and it also reaches every other + # consumer of `cell_size()`. Not fixed here because `_get_mesh_sizes` also + # feeds `get_min_radius`, the adaptivity metrics and the free-surface + # relaxation, and it needs its own benchmarking. + # Guard/measurement: tests/parallel/test_1069_boundary_normal_parallel.py + # (_nitsche_annulus_diagnostics docstring records the numbers). radii = numpy.asarray(self._radii).reshape(-1) # Empty partition (no local cells): nothing to fill on this rank. if radii.size == 0 or var.data.shape[0] == 0: @@ -3895,17 +4063,47 @@ def _do_move(): # geometry (the JIT reads the variable's .data, which would otherwise # hold the setup-time normal). Re-assemble each cached boundary normal. if getattr(self, "_boundary_normal_vars", None): + _bn_comm = self.dm.comm.tompi4py() for _nm, _var in list(self._boundary_normal_vars.items()): + # The outcome is decided COLLECTIVELY, not per rank. The callee + # already all-reduces its own failure flag and raises on every rank + # or none, so this is belt and braces for anything raised OUTSIDE + # its guarded region — but it is what makes "every rank takes the + # same branch" a property of this loop rather than an inherited + # assumption. The all-reduce is reached on the exception path too; + # that is the whole point. + _bn_failed, _bn_exc = 0, None try: self._assemble_boundary_normal(_var, _nm) - except Exception: + except Exception as _e: + _bn_failed, _bn_exc = 1, _e + if _bn_comm.size > 1: + _bn_failed = _bn_comm.allreduce(_bn_failed) + if _bn_failed: + _exc = _bn_exc if _bn_exc is not None else RuntimeError( + "failed on another rank") # Sanctioned swallow: a normal refresh can fail for a # boundary whose label vanished from the current DM # (e.g. after region extraction); the deform itself is # complete and must not be rolled back for one BC aid. # Consequence of skipping: that Nitsche/penalty BC - # keeps its setup-time normal until next refresh. - pass + # keeps its setup-time normal until next refresh — which is + # why it is a WARNING and not silence. A silent skip here + # leaves a stale constraint direction on a moved boundary. + # + # SAFE TO SWALLOW because _assemble_boundary_normal is + # COLLECTIVE (it completes the facet sum across ranks, #564) + # and its failures are collective too: it all-reduces its own + # failure flag and raises on EVERY rank or none. So what + # arrives here is already symmetric and every rank takes this + # branch together. Do not weaken that contract — a rank-local + # raise out of that routine hangs the job here rather than + # failing it. The dict is built by a collective accessor, so + # every rank iterates the same boundaries in the same order. + uw.mpi.pprint( + f"[mesh.deform] WARNING: could not refresh the boundary " + f"normal for {_nm!r} ({type(_exc).__name__}: {_exc}); BCs " + f"that captured it keep their pre-deform direction.") # Likewise refresh the local cell-size field (Nitsche penalty scaling) # so its cell-constant data tracks the deformed geometry. if getattr(self, "_cell_size_variable", None) is not None: diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index aad5153f0..13824a9d6 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/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) diff --git a/src/underworld3/utilities/facet_normals.py b/src/underworld3/utilities/facet_normals.py new file mode 100644 index 000000000..057543e6b --- /dev/null +++ b/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 diff --git a/src/underworld3/utilities/fault_contact.py b/src/underworld3/utilities/fault_contact.py index 1d3e54f8d..951b4ad2f 100644 --- a/src/underworld3/utilities/fault_contact.py +++ b/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 diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 47dfed766..640e0864d 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/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): diff --git a/tests/parallel/serial_reference.py b/tests/parallel/serial_reference.py new file mode 100644 index 000000000..ae147c3a2 --- /dev/null +++ b/tests/parallel/serial_reference.py @@ -0,0 +1,213 @@ +"""Compute a parallel test's OWN diagnostic at np=1, in THIS environment, so the +partition-independence assertion compares two runs of the same code on the same host. + +Why this exists +--------------- +The partition tests used to compare against a constant recorded on a developer's +machine, and their failure messages said "differs serial vs np=N" when what they +actually measured was "differs from a number recorded elsewhere". Those are not the +same statement, and conflating them cost a day of investigation: seven assertions +failed on CI, were read as a partition-dependence family, and four of them turned out +to be the mesh. Running ``test_1064``'s own annulus diagnostic at both rank counts on +ONE CI host gives np=1 ``1.897329151623790e-02`` and np=2 ``1.897329151623740e-02`` — +agreement to the 13th significant figure — while BOTH differ from the recorded golden +by the same +1.676e-04, because gmsh builds a different triangulation on the Linux +runner and the tests are ``mpi(min_size=2)`` so CI never ran np=1 to notice. + +A self-referential comparison is immune to that: whatever mesh the host generates, both +sides of the comparison use it. + +How +--- +Rank 0 spawns a plain single-rank Python running the test module's own ``__main__``, +which prints one ``SERIALREF `` line; the payload is broadcast to every rank. The +child's environment is scrubbed of the launcher's MPI variables — inherited ``OMPI_*`` / +``PMIX_*`` make the child believe it is a member of the parent's job, and it then hangs +or aborts instead of running as a singleton. + +The mesh is read from the same gmsh cache the parent uses, so parent and child are the +same triangulation by construction; the fingerprint is carried through anyway and +reported on failure, so the day that stops being true it says so. + +COLLECTIVE — every rank must call :func:`serial_reference` (rank 0 runs the child, the +others wait in the broadcast). +""" +import json +import os +import subprocess +import sys + +import numpy as np + +import underworld3 as uw + +# Launcher variables that would make a spawned singleton try to join the parent job. +_MPI_ENV_PREFIXES = ("OMPI_", "PMIX_", "PMI_", "MPI_", "HYDRA_", "I_MPI_", "SLURM_") + +_CACHE = {} + + +def mesh_fingerprint(mesh): + """A partition-independent identity for the mesh a diagnostic ran on: ``(global + cell count, ∫1 dV)``. A different triangulation moves both (the discretised volume + of a curved domain moves ~0.1 % between triangulations); a different partition of + the same triangulation moves neither. + + OWNED cells only. The overlap layer puts the same cell on more than one rank, so a + plain sum of the local cell-stratum sizes grows with the rank count and would make + the fingerprint report a partition as if it were a mesh change.""" + import sympy + + dm = mesh.dm + ghosts = set() + if uw.mpi.size > 1: + _nroots, local, _remote = dm.getPointSF().getGraph() + if local is not None: + ghosts = set(int(point) for point in local) + cell_start, cell_end = dm.getHeightStratum(0) + owned = sum(1 for c in range(cell_start, cell_end) if c not in ghosts) + cells = int(uw.mpi.comm.allreduce(owned)) + volume = float(uw.maths.Integral(mesh, sympy.Integer(1)).evaluate()) + return [cells, volume] + + +def serial_reference(module_file, kind, timeout=600): + """Run ``python `` as a single-rank child and return the JSON + payload it printed on its ``SERIALREF`` line. Cached per (module, kind) within the + process. Raises with the child's output if it did not produce one. + + COLLECTIVE. Rank 0 runs the child, everyone waits in the broadcast — so rank 0 + must reach that broadcast on EVERY path. ``_run_child`` therefore catches + everything and returns the failure as a string rather than raising: an + ``OSError`` from ``subprocess.run``, a truncated ``SERIALREF`` line, a + ``MemoryError`` on ``capture_output`` would otherwise leave rank 0 unwinding out + of here while every other rank sits in ``MPI_Bcast`` (busy-polling, a core each) + until pytest-timeout fires. + + The default timeout is deliberately SHORTER than the ``pytest.mark.timeout`` the + test files carry (600 s vs 900 s), so a stuck child is reported as a stuck child + instead of being overtaken by the outer timeout and reported as a stuck test. + """ + key = (os.path.abspath(module_file), kind) + if key in _CACHE: + return _CACHE[key] + + payload = _run_child(module_file, kind, timeout) if uw.mpi.rank == 0 else None + payload = uw.mpi.comm.bcast(payload, root=0) + if isinstance(payload, str): + raise RuntimeError(payload) + _CACHE[key] = payload + return payload + + +def _run_child(module_file, kind, timeout): + """Never raises — see :func:`serial_reference`. Returns the payload dict, or a + string describing what went wrong.""" + name = os.path.basename(module_file) + try: + env = {k: v for k, v in os.environ.items() + if not k.startswith(_MPI_ENV_PREFIXES)} + proc = subprocess.run( + [sys.executable, "-u", os.path.abspath(module_file), kind], + env=env, capture_output=True, text=True, timeout=timeout) + except subprocess.TimeoutExpired: + return f"serial reference for {name}:{kind} timed out after {timeout}s" + except Exception as exc: # noqa: BLE001 - see the docstring + return (f"serial reference for {name}:{kind} could not be launched — " + f"{type(exc).__name__}: {exc}") + try: + for line in proc.stdout.splitlines(): + if line.startswith("SERIALREF "): + return json.loads(line[len("SERIALREF "):]) + except Exception as exc: # noqa: BLE001 - see the docstring + return (f"serial reference for {name}:{kind} printed an unreadable " + f"SERIALREF line — {type(exc).__name__}: {exc}") + return (f"serial reference for {name}:{kind} printed no " + f"SERIALREF line (rc={proc.returncode})\n" + f"--- stdout tail ---\n{proc.stdout[-2000:]}\n" + f"--- stderr tail ---\n{proc.stderr[-2000:]}") + + +def emit(values, fingerprint): + """Print the ``SERIALREF`` line a test module's ``__main__`` owes its parallel + twin. Rank-0 only, so it is safe to call unconditionally.""" + if uw.mpi.rank == 0: + print("SERIALREF " + json.dumps( + {"values": [float(v) for v in np.atleast_1d(values)], + "fingerprint": [float(f) for f in fingerprint]})) + + +def _same_mesh(left, right): + """Do two fingerprints describe the same triangulation? Cell count exactly, volume + to 1e-12 relative (it is a sum of the same element volumes in a partition-dependent + ORDER, so the last couple of bits move; a different triangulation moves it by ~1e-3 + relative, nine orders away).""" + return (int(left[0]) == int(right[0]) + and np.isclose(left[1], right[1], rtol=1e-12, atol=0)) + + +def _fp(fingerprint): + return f"cells={fingerprint[0]:.0f} vol={fingerprint[1]:.12g}" + + +def compare(values, reference, rtols, labels, fingerprint, what): + """Assert each of ``values`` matches the serial reference within its ``rtols``, + and say what moved — including both mesh fingerprints, so a host/mesh difference + reads as a mesh difference instead of as a physics regression. + + The fingerprints are ASSERTED equal, not merely reported. The np=1 child reads the + same gmsh cache as the parallel parent, so a mismatch means something has broken + that assumption (a concurrent run regenerating the cache, a fingerprint that is not + partition-independent after all) and every number below it would be meaningless. + """ + ref_values = reference["values"] + ref_fp = reference["fingerprint"] + assert len(values) == len(ref_values), ( + f"{what}: serial reference has {len(ref_values)} values, this run produced " + f"{len(values)}") + assert _same_mesh(fingerprint, ref_fp), ( + f"{what}: the np=1 reference ran on a DIFFERENT mesh from this np={uw.mpi.size} " + f"run — np=1 [{_fp(ref_fp)}] vs np={uw.mpi.size} [{_fp(fingerprint)}]. The " + f"comparison below would be measuring the mesh, not the partition.") + fp_note = f" [mesh {_fp(fingerprint)}]" + for value, ref, rtol, label in zip(values, ref_values, rtols, labels): + assert np.isclose(value, ref, rtol=rtol, atol=0), ( + f"{what}: {label} is partition dependent — np=1 {ref!r} vs " + f"np={uw.mpi.size} {value!r} (rtol {rtol:g}){fp_note}") + + +def accuracy_anchor(values, anchor, fingerprint, labels, what, rtol=1e-2): + """Assert the ABSOLUTE answer against a recorded constant, GATED on the mesh. + + Partition independence and accuracy are two different claims and the self-referential + comparison above only makes the first. A rotated constraint that stopped constraining + equally on every rank, an FMG hierarchy that converged to the wrong answer, a + physics benchmark coefficient that drifted — all of those pass ``compare`` and are + caught only by a number recorded when the result was known good. + + The reason those constants were removed is real: they are host-specific, because + gmsh triangulates differently on different platforms, and a mismatch then reads as a + physics regression. The fix is the one the #564 investigation actually recommended — + keep the constant, and put the mesh fingerprint in front of it. On the host the + anchor was recorded on this is a live accuracy gate; on any other mesh it SKIPS, + loudly, instead of failing for the wrong reason. + + ``rtol`` is deliberately loose (1 %). This is not a reproducibility check — that is + ``compare``'s job, three to eight orders tighter. This one only has to notice that + the answer has become a different answer. + """ + import pytest + + if not _same_mesh(fingerprint, anchor["fingerprint"]): + pytest.skip( + f"{what}: accuracy anchor was recorded on a different mesh " + f"[{_fp(anchor['fingerprint'])}] from this host's " + f"[{_fp(fingerprint)}] — gmsh triangulates differently across " + f"platforms. Partition independence is still asserted; only the " + f"absolute value is skipped.") + for value, ref, label in zip(values, anchor["values"], labels): + assert np.isclose(value, ref, rtol=rtol, atol=0), ( + f"{what}: {label} has MOVED from its recorded value on the same mesh " + f"[{_fp(fingerprint)}] — {ref!r} recorded, {value!r} now " + f"(rtol {rtol:g}). This is an accuracy regression, not a partition " + f"effect: the mesh is identical and the partition check passed.") diff --git a/tests/parallel/test_1063_constrained_freeslip_parallel.py b/tests/parallel/test_1063_constrained_freeslip_parallel.py index 843f99a08..260ee4684 100644 --- a/tests/parallel/test_1063_constrained_freeslip_parallel.py +++ b/tests/parallel/test_1063_constrained_freeslip_parallel.py @@ -5,18 +5,33 @@ surgery, so the GLOBAL system — and hence the velocity solve and the gauge-invariant boundary traction — are partition-independent. This test verifies that, for both isotropic and transverse-isotropic rheology, the -parallel solve reproduces the serial reference to a **tight tolerance** (the +parallel solve reproduces its OWN np=1 answer to a **tight tolerance** (the residual difference is the parallel reduction order, not the solver): - * the velocity L2 norm (``∫ v·v``), and - * the MEAN-STRIPPED boundary topography (``∫(h - h̄)²`` on the constrained - boundary), computed via ``solver.topography(boundary, reference="mean")``. + * the velocity L2 norm (``int v.v``), and + * the MEAN-STRIPPED boundary topography (``int (h - hbar)^2`` on the + constrained boundary), via ``solver.topography(boundary, reference="mean")``. -The raw multiplier ``h`` carries the ``[p,λ]`` gauge constant — the solver lands -on a partition-dependent representative of it — so only the mean-stripped +The raw multiplier ``h`` carries the ``[p,lambda]`` gauge constant — the solver +lands on a partition-dependent representative of it — so only the mean-stripped (physical) topography is compared. All diagnostics use the parallel-safe ``uw.maths.Integral`` / ``BdIntegral`` reductions (no direct mpi4py). +The np=1 side is computed by running this file's own ``__main__`` as a single-rank +child in the same environment (``serial_reference``), NOT read from a constant +recorded on a developer's machine. That distinction matters twice over here: a +self-referential comparison is the property these tests claim to be testing, and +comparing against a stored constant is what made #564 look like seven instances of +one defect when four of them were the host's gmsh building a different mesh. + +What this used to measure, and no longer does. With a rank-local boundary normal +(#564, fixed) the ``[iso]`` velocity L2 was 6.194547793939e-01 at np=1 against +5.982807168537e-01 at np=2 — **3.4 %**, #495's number to every digit — with the +topography 2.4 % adrift; ``[ti]`` moved 0.34 % at np=2 and 2.8 % at np=4. The +default constraint direction is ``mesh.boundary_normal(boundary)``, which was +assembled from each rank's own facets. It now agrees to 1.2e-10 (iso) and 4.1e-10 +(ti) across np=1,2,3,4. + Run with: mpirun -n 2 python -m pytest --with-mpi \\ tests/parallel/test_1063_constrained_freeslip_parallel.py @@ -30,20 +45,48 @@ import underworld3 as uw -pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(180)] - -# SERIAL (np=1) reference diagnostics — the partition-independent ground truth. -# Computed once with `python {iso,ti}`; the parallel run must match. -# (velocity L2, mean-stripped boundary topography) for the annulus problem below. -GOLDEN = { - "iso": (6.194547793955e-01, 3.786068778041e+01), - "ti": (3.925981604039e-01, 3.707799837159e+01), +from serial_reference import ( + accuracy_anchor, compare, emit, mesh_fingerprint, serial_reference) + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(600)] + +# ABSOLUTE accuracy anchors, gated on the mesh fingerprint. +# +# The self-referential comparison below proves the answer does not depend on the +# PARTITION. It says nothing about whether the answer is RIGHT — a solve that broke +# identically at every rank count would sail through it. These are the pre-#568 +# goldens, kept for that second job and gated so that a host whose gmsh triangulates +# differently SKIPS the accuracy claim rather than failing it. That is what the #564 +# investigation recommended and it costs nothing. +# +# rtol is 1e-2 on purpose: this is not a reproducibility gate (that is `compare`, six +# to eight orders tighter), it only has to notice that the answer became a different +# answer. Recompute the fingerprint with `python `. +_ANCHOR_MESH = [434, 2.355685412534114] # Annulus(1.0, 0.5, cellSize=0.12) +ANCHORS = { + "iso": {"fingerprint": _ANCHOR_MESH, + "values": (6.194547793955e-01, 3.786068778041e+01)}, + "ti": {"fingerprint": _ANCHOR_MESH, + "values": (3.925981604039e-01, 3.707799837159e+01)}, } +# (velocity L2, mean-stripped topography). The raw mean pressure is NOT anchored: it is +# pinned to ~0, and a relative gate on a number near machine zero measures nothing — +# it has its own absolute assertion in the test. +ANCHOR_GAUGE = {"fingerprint": _ANCHOR_MESH, + "values": (6.194547487092e-01, 3.781793823254e+01)} + +# Velocity reproduces to the parallel reduction order. The topography is read off +# the multiplier, whose [p,lambda] Schur sub-block grinds into its 200-iteration cap +# on this problem (a converged SNES over a capped inner block — see #564's +# investigation), so it reproduces two to three orders less tightly. Both gates sit +# at least three orders below the pre-fix move they exist to catch (3.4 % and 2.4 %). +_RTOL_VELOCITY = 1.0e-8 +_RTOL_TOPOGRAPHY = 1.0e-5 def _solve_diagnostics(kind): - """Build + solve the constrained free-slip annulus; return partition- - independent diagnostics (L2 velocity, mean-stripped boundary topography).""" + """Build + solve the constrained free-slip annulus; return partition-independent + diagnostics ``((L2 velocity, mean-stripped boundary topography), fingerprint)``.""" mesh = uw.meshing.Annulus(radiusOuter=1.0, radiusInner=0.5, cellSize=0.12, qdegree=4) v = uw.discretisation.MeshVariable("U", mesh, 2, degree=2) @@ -63,6 +106,8 @@ def _solve_diagnostics(kind): st.constitutive_model.Parameters.shear_viscosity_0 = 1.0 st.tolerance = 1.0e-8 st.add_essential_bc((0.0, 0.0), "Lower") + # NO normal= : this is the path that reads mesh.boundary_normal(), which is what + # #564 was about. Passing an analytic normal here would test something else. h = st.add_constraint_bc(0.0, "Upper") st.bodyforce = 1.0e2 * sympy.sin(3 * sympy.atan2(X[1], X[0])) * unit_r st.solve(zero_init_guess=True) @@ -73,24 +118,22 @@ def _solve_diagnostics(kind): topo_fn = st.topography("Upper", reference="mean") topo = float(np.sqrt(uw.maths.BdIntegral( mesh=mesh, fn=topo_fn ** 2, boundary="Upper").evaluate())) - return L2, topo - - -# Serial (np=1) reference for the gauge-reproducibility test below. The enclosed -# iso problem with an ACTIVE pressure null space exercises the automatic pressure -# gauge (auto_pressure_gauge, default on): the constant pressure and constant -# multiplier are both gauge-free, so without a pin the solver lands on a -# partition-dependent level for each. The auto gauge pins the raw PRESSURE -# reproducibly (the raw multiplier keeps its own gauge freedom — topography is -# read gauge-invariantly via reference="mean"). Velocity is physics-neutral. -# (velocity L2, raw meanP, MEAN-STRIPPED boundary topography). -# Recompute with `python gauge`. -GOLDEN_GAUGE = (6.194547487092e-01, -3.847657609797e-10, 3.781793823254e+01) + return (L2, topo), mesh_fingerprint(mesh) def _solve_gauge_diagnostics(): """Enclosed constrained annulus with an active pressure null space (so the - automatic pressure gauge fires); return the gauge-relevant diagnostics.""" + automatic pressure gauge fires); return the gauge-relevant diagnostics + ``((L2 velocity, raw mean pressure, mean-stripped topography), fingerprint)``. + + The enclosed iso problem with an ACTIVE pressure null space exercises the + automatic pressure gauge (``auto_pressure_gauge``, default on): the constant + pressure and the constant multiplier are both gauge-free, so without a pin the + solver lands on a partition-dependent level for each. The auto gauge pins the raw + PRESSURE reproducibly; the raw multiplier keeps its own gauge freedom, which is + why topography is read gauge-invariantly via ``reference="mean"``. Velocity is + physics-neutral under the gauge. + """ mesh = uw.meshing.Annulus(radiusOuter=1.0, radiusInner=0.5, cellSize=0.12, qdegree=4) v = uw.discretisation.MeshVariable("Ug", mesh, 2, degree=2) @@ -118,18 +161,9 @@ def _solve_gauge_diagnostics(): topo_fn = st.topography("Upper", reference="mean") topoL2 = float(np.sqrt(uw.maths.BdIntegral( mesh=mesh, fn=topo_fn ** 2, boundary="Upper").evaluate())) - return L2, meanP, topoL2 - - -@pytest.mark.xfail( - reason="#564: free-slip solves are partition dependent. CI measures velocity L2 " - "0.6194547487092 serial vs 0.6194402844556 at np=2 (2.3e-05). PRE-EXISTING " - "and not caused by #560/#561: the same seven assertions fail with numbers " - "identical to every digit at #561's merge base with only the " - "scripts/test.sh test_10*py line enabled, which is how they became visible " - "at all — this whole batch had never run in CI. Passes locally on " - "macOS/arm64, so strict=False; see #564 for the full table.", - strict=False) + return (L2, meanP, topoL2), mesh_fingerprint(mesh) + + def test_constrained_raw_gauge_partition_independent(): """With the automatic pressure gauge on (default), the RAW mean pressure is partition-independent (pinned to ~0), the velocity stays bit-identical (the @@ -137,58 +171,46 @@ def test_constrained_raw_gauge_partition_independent(): The raw multiplier level is NOT asserted reproducible — it keeps an independent gauge freedom the pressure pin does not touch (use reference="mean").""" - L2, meanP, topo = _solve_gauge_diagnostics() - L2_ref, meanP_ref, topo_ref = GOLDEN_GAUGE - # Velocity is physics-neutral under the gauge; the residual ~1e-9 spread is - # parallel round-off in the pressure-null-space projection (this enclosed - # config carries the null space, unlike the half-pinned case above), well - # below any real partition effect. - assert np.isclose(L2, L2_ref, rtol=1e-8, atol=0), ( - f"velocity L2 differs serial vs np={uw.mpi.size}: {L2_ref} vs {L2}") - # meanP is pinned to ~0 on the gauge boundary; compare on an absolute scale - # (a relative tolerance is meaningless near machine zero). - assert np.isclose(meanP, meanP_ref, rtol=0, atol=1e-6), ( - f"raw mean pressure differs serial vs np={uw.mpi.size}: " - f"{meanP_ref} vs {meanP}") - assert np.isclose(topo, topo_ref, rtol=1e-6, atol=0), ( - f"mean-stripped topography differs serial vs np={uw.mpi.size}: " - f"{topo_ref} vs {topo}") - - -@pytest.mark.xfail( - reason="#564 (of which #495 is one member): free-slip solves are partition " - "dependent. CI measures velocity L2 [iso] 0.6194547793955 serial vs " - "0.6107410846031 at np=2 (1.4%) and [ti] 0.3925981604039 vs " - "0.3937854671587 (0.3%), against the 1e-9 this asserts. PRE-EXISTING and " - "not caused by #560/#561: the same numbers reproduce to every digit at " - "#561's merge base with only the scripts/test.sh test_10*py line enabled, " - "and Stokes_Constrained never calls _boundary_velocity_nodes. #564 records " - "that the ROTATED path is affected too, so this is a family rather than a " - "single solver's bug. Passes locally on macOS/arm64, so strict=False.", - strict=False) + (L2, meanP, topo), fingerprint = _solve_gauge_diagnostics() + reference = serial_reference(__file__, "gauge") + # meanP is pinned to ~0 on the gauge boundary, so the claim is an ABSOLUTE one + # and needs no reference at all — a relative comparison of two numbers near + # machine zero measures nothing. The other two go against the np=1 run. + assert abs(meanP) < 1e-6, ( + f"raw mean pressure is not pinned at np={uw.mpi.size}: {meanP!r}") + labels = ("velocity L2", "mean-stripped topography") + compare((L2, topo), {"values": [reference["values"][0], reference["values"][2]], + "fingerprint": reference["fingerprint"]}, + rtols=(_RTOL_VELOCITY, _RTOL_TOPOGRAPHY), labels=labels, + fingerprint=fingerprint, what="constrained raw-gauge annulus") + accuracy_anchor((L2, topo), ANCHOR_GAUGE, fingerprint, labels, + what="constrained raw-gauge annulus") + + @pytest.mark.parametrize("kind", ["iso", "ti"]) def test_constrained_freeslip_partition_independent(kind): - """The parallel solve must reproduce the serial reference: velocity bit- - identical, and the gauge-fixed (mean-stripped) topography bit-identical.""" - L2_par, topo_par = _solve_diagnostics(kind) - L2_ref, topo_ref = GOLDEN[kind] - assert np.isclose(L2_par, L2_ref, rtol=1e-9, atol=0), ( - f"[{kind}] velocity L2 differs serial vs np={uw.mpi.size}: " - f"{L2_ref} vs {L2_par}") - assert np.isclose(topo_par, topo_ref, rtol=1e-6, atol=0), ( - f"[{kind}] mean-stripped topography differs serial vs np={uw.mpi.size}: " - f"{topo_ref} vs {topo_par}") + """The parallel solve must reproduce its own np=1 answer: velocity to the + parallel reduction order, and the gauge-fixed (mean-stripped) topography to the + multiplier solve's reproducibility.""" + values, fingerprint = _solve_diagnostics(kind) + labels = ("velocity L2", "mean-stripped topography") + compare(values, serial_reference(__file__, kind), + rtols=(_RTOL_VELOCITY, _RTOL_TOPOGRAPHY), labels=labels, + fingerprint=fingerprint, what=f"constrained free-slip [{kind}]") + accuracy_anchor(values, ANCHORS[kind], fingerprint, labels, + what=f"constrained free-slip [{kind}]") if __name__ == "__main__": - # Recompute the serial GOLDEN reference: `python {iso,ti,gauge}`. + # Single-rank child of the parallel run (see serial_reference), and a + # human-readable recompute: `python {iso,ti,gauge}`. import sys _kind = sys.argv[1] if len(sys.argv) > 1 else "iso" if _kind == "gauge": - _L2, _meanP, _topo = _solve_gauge_diagnostics() - if uw.mpi.rank == 0: - print(f"DIAG_GAUGE {_L2:.12e} {_meanP:.12e} {_topo:.12e}") + _values, _fingerprint = _solve_gauge_diagnostics() + emit(_values, _fingerprint) + uw.mpi.pprint("DIAG_GAUGE " + " ".join(f"{v:.12e}" for v in _values)) else: - _L2, _topo = _solve_diagnostics(_kind) - if uw.mpi.rank == 0: - print(f"DIAG {_L2:.12e} {_topo:.12e}") + _values, _fingerprint = _solve_diagnostics(_kind) + emit(_values, _fingerprint) + uw.mpi.pprint("DIAG " + " ".join(f"{v:.12e}" for v in _values)) diff --git a/tests/parallel/test_1064_rotated_freeslip_parallel.py b/tests/parallel/test_1064_rotated_freeslip_parallel.py index 63a935a11..aeed50d4d 100644 --- a/tests/parallel/test_1064_rotated_freeslip_parallel.py +++ b/tests/parallel/test_1064_rotated_freeslip_parallel.py @@ -8,7 +8,7 @@ With that fixed, the whole global system (and hence the velocity solve and the wall-normal leakage) is partition-independent. -This test verifies that the parallel solve reproduces the serial reference to a tight +This test verifies that the parallel solve reproduces ITS OWN np=1 answer to a tight tolerance for two geometries: * **box** — 4-wall rotated free-slip on axis-aligned walls (GAMG velocity block); the @@ -22,6 +22,26 @@ (no rank-local ``v.data``, which is per-partition, and no ``uw.function.evaluate`` on arbitrary points, which deadlocks np>1). +Every "must match serial" assertion below compares against a np=1 run of THIS FILE, +computed in THIS environment by ``serial_reference`` — not against a constant recorded +on a developer's machine. The distinction is not cosmetic. Five of these tests carried +hardcoded goldens and messages reading "differs serial vs np=N", and on CI they failed: +annulus velocity L2 recorded 1.897011154231e-02, measured 1.897329151624e-02. Running +the SAME diagnostic at BOTH rank counts on one CI host gives np=1 1.897329151623790e-02 +and np=2 1.897329151623740e-02 — agreement to the 13th significant figure, with the +leakages identical to every digit — and both differ from the recorded golden by exactly +the same +1.676e-04. The rotated solve was never partition dependent. gmsh builds a +different triangulation on the Linux runner than on macOS/arm64, the goldens were +recorded on macOS, and because these tests are ``mpi(min_size=2)`` CI never ran np=1 to +notice. A day of investigation went into a defect that did not exist, and the +misleading part was the assertion message. + +Where a genuine ABSOLUTE check is wanted — the analytic SolCx velocity error, the +sigma_nn accuracy against the exact solution, the iteration-count ceiling — it is kept +and labelled as an accuracy check, with a tolerance loose enough to survive a +cross-host mesh change. Those are NOT partition-independence checks and must not be +conflated with them again. + Run with: mpirun -n 2 python -m pytest --with-mpi \\ tests/parallel/test_1064_rotated_freeslip_parallel.py @@ -37,41 +57,53 @@ from underworld3.function import analytic as A from underworld3.utilities import custom_mg -pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(180)] - -# SERIAL (np=1) reference diagnostics — the partition-independent ground truth. -# Recompute with `python {box,annulus}`. -# box: (velocity L2, analytic velocity error) -# annulus: (velocity L2, radial-leakage L2 on Lower arc, radial-leakage L2 on Upper arc) -GOLDEN_BOX = (1.275109036912e-03, 1.529545e-05) -GOLDEN_ANNULUS = (1.897011154231e-02, 4.563841e-05, 9.341699e-06) -# annulus driven by CUSTOM GEOMETRIC FMG on the velocity block (nested hierarchy): -# (velocity L2, radial-leakage L2 on Lower arc, radial-leakage L2 on Upper arc) -GOLDEN_ANNULUS_FMG = (1.906961759626e-02, 5.428193e-06, 1.177002e-06) -# 3D spherical shell (free-slip both boundaries, all 3 rotation nullspace modes): -# velocity L2. Recompute with `python spherical3d`. -GOLDEN_SPHERICAL3D = 4.069689334228e-03 -# Zhong l=2 topography coefficients recovered from the 3D rotated-constraint -# reaction: surface (all, vertices, midpoints), CMB (all, vertices, midpoints). -# Recompute with `python spherical3d_topo`. -GOLDEN_SPHERICAL3D_TOPO = ( - 4.149689252074e-01, - 3.952301937705e-01, - 4.215939953379e-01, - 7.932177563075e-01, - 8.426041179682e-01, - 7.762363224500e-01, -) -# NONLINEAR (power-law) box with rotated free-slip through the manual Newton loop -# (consistent tangent): (velocity L2, nonlinear iteration count — the number of -# Newton increments solved, == len(ksp_its); this solve exits on the step-norm -# test after its 8th increment). Recompute `python nonlinear`. -GOLDEN_BOX_NONLINEAR = (8.069396188270e-04, 8) -# box sigma_nn (boundary_normal_traction on Top, default lumped mass) vs analytic SolCx -# sigma_yy, whole boundary: (relL2, |corr|). Recompute with `python sigma`. -GOLDEN_BOX_SIGMA = (5.554578e-02, 0.998466) -# dynamic_topography field: BdIntegral L2 of h over Top. Recompute `python topo`. -GOLDEN_TOPO_BDL2 = 2.553916470e-01 +from serial_reference import ( + accuracy_anchor, compare, emit, mesh_fingerprint, serial_reference) + +# The timeout covers the np=1 child each partition test spawns as well as the +# parallel solve itself. +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(900)] + +# ABSOLUTE accuracy anchors, gated on the mesh fingerprint. +# +# `compare` proves the answer does not depend on the PARTITION; it says nothing about +# whether the answer is RIGHT. A rotated constraint that stopped constraining equally +# on every rank, an FMG hierarchy converging to the wrong place, a Zhong l=2 benchmark +# coefficient drifting — all of those pass a self-referential test. These are the +# pre-#568 goldens, kept for that second job, and gated so a host whose gmsh +# triangulates differently SKIPS the accuracy claim instead of failing it (which is +# what made them misleading in the first place — see the module docstring). +# +# rtol is 1e-2: not a reproducibility gate, just "is this still the same answer". +# Fingerprints are (owned cell count, integral 1 dV); recompute any line with +# `python `. +_MESH_BOX24 = [576, 1.0000000000000004] # StructuredQuadBox(24, 24) +_MESH_BOX8 = [64, 0.9999999999999999] # StructuredQuadBox(8, 8) +_MESH_ANNULUS = [600, 2.356187202481425] # Annulus(0.5, 1.0, cellSize=0.1) +_MESH_ANNULUS_FMG = [2304, 2.356078287527854] # cellSize=0.2, refined twice +_MESH_SHELL = [1417, 3.4521585981151097] # SphericalShell(0.55, 1.0, cs=0.25) +_MESH_SHELL_INT = [2016, 3.4521585981151106] # SphericalShellInternalBoundary + +ANCHORS = { + "box": {"fingerprint": _MESH_BOX24, + "values": (1.275109036912e-03,)}, + "annulus": {"fingerprint": _MESH_ANNULUS, + "values": (1.897011154231e-02, 4.563841e-05, 9.341699e-06)}, + "annulus_fmg": {"fingerprint": _MESH_ANNULUS_FMG, + "values": (1.906961759626e-02, 5.428193e-06, 1.177002e-06)}, + "spherical3d": {"fingerprint": _MESH_SHELL, + "values": (4.069689334228e-03,)}, + "spherical3d_topo": {"fingerprint": _MESH_SHELL_INT, + "values": (4.149689252074e-01, 3.952301937705e-01, + 4.215939953379e-01, 7.932177563075e-01, + 8.426041179682e-01, 7.762363224500e-01)}, + "nonlinear": {"fingerprint": _MESH_BOX8, + "values": (8.069396188270e-04, 8)}, + "sigma": {"fingerprint": _MESH_BOX24, + "values": (5.554578e-02, 0.998466)}, + "topo": {"fingerprint": _MESH_BOX24, + "values": (2.553916470e-01,)}, +} def _wrap(dm, m0): @@ -105,7 +137,7 @@ def _box_diagnostics(): L2 = float(np.sqrt(uw.maths.Integral(mesh, v.sym.dot(v.sym)).evaluate())) verr = float(sol.velocity_error(v)) - return L2, verr + return (L2, verr), mesh_fingerprint(mesh) def _annulus_diagnostics(): @@ -137,7 +169,7 @@ def _annulus_diagnostics(): mesh=mesh, fn=vr**2, boundary="Lower").evaluate())) leak_up = float(np.sqrt(uw.maths.BdIntegral( mesh=mesh, fn=vr**2, boundary="Upper").evaluate())) - return L2, leak_lo, leak_up + return (L2, leak_lo, leak_up), mesh_fingerprint(mesh) def _spherical3d_diagnostics(): @@ -167,7 +199,8 @@ def _spherical3d_diagnostics(): L2 = float(np.sqrt(uw.maths.Integral(mesh, v.sym.dot(v.sym)).evaluate())) info = s._rotated_freeslip_info # the unified rotated loop reports one KSP count per Newton increment - return L2, max(info["ksp_its"]), int(info["ksp_reason"]) + return ((L2,), mesh_fingerprint(mesh), + max(info["ksp_its"]), int(info["ksp_reason"])) def _spherical3d_topography_diagnostics(cell_size=0.25): @@ -233,7 +266,7 @@ def fit(mask): return ( *harmonic_coefficients("Upper", 1.0), *harmonic_coefficients("Lower", -1.0), - ) + ), mesh_fingerprint(mesh) def _annulus_fmg_diagnostics(mesh_owned=False): @@ -285,7 +318,7 @@ def _annulus_fmg_diagnostics(mesh_owned=False): mesh=fine, fn=vr**2, boundary="Lower").evaluate())) leak_up = float(np.sqrt(uw.maths.BdIntegral( mesh=fine, fn=vr**2, boundary="Upper").evaluate())) - return L2, leak_lo, leak_up + return (L2, leak_lo, leak_up), mesh_fingerprint(fine) def _box_nonlinear_diagnostics(): @@ -316,7 +349,8 @@ def _box_nonlinear_diagnostics(): s.solve() L2 = float(np.sqrt(uw.maths.Integral(mesh, v.sym.dot(v.sym)).evaluate())) - return L2, int(s._rotated_freeslip_info["nonlinear_iterations"]) + return ((L2, int(s._rotated_freeslip_info["nonlinear_iterations"])), + mesh_fingerprint(mesh)) def _box_sigma_diagnostics(): @@ -361,7 +395,7 @@ def _box_sigma_diagnostics(): S = S if corr >= 0 else -S relL2 = float(np.linalg.norm(S - syy) / np.linalg.norm(syy)) result = (relL2, abs(corr), len(X)) - return comm.bcast(result, root=0) + return comm.bcast(result, root=0), mesh_fingerprint(mesh) def _box_topography_bdl2(): @@ -386,65 +420,59 @@ def _box_topography_bdl2(): s.petsc_options["snes_type"] = "ksponly" s.solve() s.dynamic_topography("Top", hf, buoyancy_scale=1.0) - return float(np.sqrt(uw.maths.BdIntegral( + bdl2 = float(np.sqrt(uw.maths.BdIntegral( mesh=mesh, fn=hf.sym[0] ** 2, boundary="Top").evaluate())) + return (bdl2,), mesh_fingerprint(mesh) def test_rotated_freeslip_box_partition_independent(): - """Box: the parallel rotated free-slip solve reproduces the serial velocity L2 and - keeps the analytic velocity error small.""" - L2, verr = _box_diagnostics() - L2_ref, verr_ref = GOLDEN_BOX - assert np.isclose(L2, L2_ref, rtol=1e-8, atol=0), ( - f"box velocity L2 differs serial vs np={uw.mpi.size}: {L2_ref} vs {L2}") - # analytic accuracy is preserved (partition may change the exact digit but not - # the order of magnitude of the SolCx error) - assert verr < 1e-3, f"box velocity error {verr:.2e} too large at np={uw.mpi.size}" - - -@pytest.mark.xfail( - reason="#564: free-slip solves are partition dependent. CI measures annulus " - "velocity L2 0.01897011154231 serial vs 0.01897329151624 at np=2 (1.7e-04). " - "PRE-EXISTING and not caused by #560/#561: the same seven assertions fail " - "with numbers identical to every digit at #561's merge base with only the " - "scripts/test.sh test_10*py line enabled, which is how they became visible " - "at all — this whole batch had never run in CI. This case passes an " - "explicit analytic normal=, and its numbers are bit-identical before and " - "after #560's nodal-normal fix, so it is definitively not that mechanism. " - "Passes locally on macOS/arm64, so strict=False; see #564 for the full " - "table.", - strict=False) + """Box: the parallel rotated free-slip solve reproduces its OWN np=1 velocity L2, + and the analytic velocity error stays small.""" + values, fingerprint = _box_diagnostics() + compare(values[:1], _reference("box", 1), rtols=(1e-8,), labels=("velocity L2",), + fingerprint=fingerprint, what="box rotated free-slip") + # ACCURACY check, not a partition check: the SolCx error is a property of the + # discretisation, and the gate is loose enough to survive a cross-host mesh. + # Asserted BEFORE the fingerprint-gated anchor: that anchor SKIPS on a host + # whose mesh differs, and this assertion must not be skipped with it. + assert values[1] < 1e-3, ( + f"box velocity error {values[1]:.2e} too large at np={uw.mpi.size}") + accuracy_anchor(values[:1], ANCHORS["box"], fingerprint, ("velocity L2",), + what="box rotated free-slip") + + def test_rotated_freeslip_annulus_partition_independent(): - """Annulus: the parallel radial free-slip solve reproduces the serial velocity L2 - and the (partition-independent) radial leakage on both arcs.""" - L2, leak_lo, leak_up = _annulus_diagnostics() - L2_ref, leak_lo_ref, leak_up_ref = GOLDEN_ANNULUS + """Annulus: the parallel radial free-slip solve reproduces its OWN np=1 velocity L2 + and radial leakage on both arcs. + + Measured on ONE CI host, this diagnostic: np=1 1.897329151623790e-02, np=2 + 1.897329151623740e-02, leakages identical to every digit. The 1.7e-04 this test + used to report was the distance to a golden recorded on a different host, not a + partition effect — see the module docstring. + """ + values, fingerprint = _annulus_diagnostics() # velocity L2 is iterative-solver-tolerance reproducible (~1e-8 rel), not the # box's 1e-10 — the annulus carries a rotation null space + gauge removal. - assert np.isclose(L2, L2_ref, rtol=1e-6, atol=0), ( - f"annulus velocity L2 differs serial vs np={uw.mpi.size}: {L2_ref} vs {L2}") - assert np.isclose(leak_lo, leak_lo_ref, rtol=1e-4, atol=0), ( - f"annulus Lower leakage differs serial vs np={uw.mpi.size}: " - f"{leak_lo_ref} vs {leak_lo}") - assert np.isclose(leak_up, leak_up_ref, rtol=1e-4, atol=0), ( - f"annulus Upper leakage differs serial vs np={uw.mpi.size}: " - f"{leak_up_ref} vs {leak_up}") + labels = ("velocity L2", "Lower radial leakage", "Upper radial leakage") + compare(values, _reference("annulus", 3), rtols=(1e-6, 1e-4, 1e-4), + labels=labels, fingerprint=fingerprint, what="annulus rotated free-slip") + # The LEAKAGE anchors are the point here: a rotated constraint that stopped + # constraining would keep every partition agreeing with every other. + accuracy_anchor(values, ANCHORS["annulus"], fingerprint, labels, + what="annulus rotated free-slip") def test_rotated_freeslip_annulus_fmg_partition_independent(): """The full stack: rotated radial free-slip on the annulus with the velocity block - driven by CUSTOM GEOMETRIC FMG (set_custom_fmg) reproduces the serial velocity L2 + driven by CUSTOM GEOMETRIC FMG (set_custom_fmg) reproduces its own np=1 velocity L2 and radial leakage in parallel — FMG x rotated x annulus x np>1.""" - L2, leak_lo, leak_up = _annulus_fmg_diagnostics() - L2_ref, leak_lo_ref, leak_up_ref = GOLDEN_ANNULUS_FMG - assert np.isclose(L2, L2_ref, rtol=1e-6, atol=0), ( - f"FMG annulus velocity L2 differs serial vs np={uw.mpi.size}: {L2_ref} vs {L2}") - assert np.isclose(leak_lo, leak_lo_ref, rtol=1e-4, atol=0), ( - f"FMG annulus Lower leakage differs serial vs np={uw.mpi.size}: " - f"{leak_lo_ref} vs {leak_lo}") - assert np.isclose(leak_up, leak_up_ref, rtol=1e-4, atol=0), ( - f"FMG annulus Upper leakage differs serial vs np={uw.mpi.size}: " - f"{leak_up_ref} vs {leak_up}") + values, fingerprint = _annulus_fmg_diagnostics() + labels = ("velocity L2", "Lower radial leakage", "Upper radial leakage") + compare(values, _reference("annulus_fmg", 3), rtols=(1e-6, 1e-4, 1e-4), + labels=labels, fingerprint=fingerprint, + what="custom-FMG annulus rotated free-slip") + accuracy_anchor(values, ANCHORS["annulus_fmg"], fingerprint, labels, + what="custom-FMG annulus rotated free-slip") def test_rotated_freeslip_mesh_owned_fmg_pickup(): @@ -456,151 +484,144 @@ def test_rotated_freeslip_mesh_owned_fmg_pickup(): The transfers are built cross-partition, so this is not implied by the serial pickup test (``tests/test_1021_mg_option_bundle.py``); the tail is attached by hand because what is under test is whether the rotated dispatch consults it, - not how ``adapt()`` produces it.""" - L2, leak_lo, leak_up = _annulus_fmg_diagnostics(mesh_owned=True) - L2_ref, leak_lo_ref, leak_up_ref = GOLDEN_ANNULUS_FMG - assert np.isclose(L2, L2_ref, rtol=1e-6, atol=0), ( - f"mesh-owned FMG annulus velocity L2 differs serial vs np={uw.mpi.size}: " - f"{L2_ref} vs {L2}") - assert np.isclose(leak_lo, leak_lo_ref, rtol=1e-4, atol=0) - assert np.isclose(leak_up, leak_up_ref, rtol=1e-4, atol=0) - - -@pytest.mark.xfail( - reason="#564: free-slip solves are partition dependent. CI measures 3-D spherical " - "velocity L2 0.004069689334228 serial vs 0.004074314572473 at np=2 " - "(1.1e-03). PRE-EXISTING and not caused by #560/#561: the same seven " - "assertions fail with numbers identical to every digit at #561's merge base " - "with only the scripts/test.sh test_10*py line enabled, which is how they " - "became visible at all — this whole batch had never run in CI. This case " - "passes an explicit analytic normal=, and its numbers are bit-identical " - "before and after #560's nodal-normal fix, so it is definitively not that " - "mechanism. Passes locally on macOS/arm64, so strict=False; see #564 for " - "the full table.", - strict=False) + not how ``adapt()`` produces it. + + The np=1 reference is the EXPLICIT-registration one: #467 is precisely the claim + that the two routes produce the same solve, so comparing the mesh-owned parallel + run against the explicitly-registered serial run asserts both properties at once. + """ + values, fingerprint = _annulus_fmg_diagnostics(mesh_owned=True) + labels = ("velocity L2", "Lower radial leakage", "Upper radial leakage") + compare(values, _reference("annulus_fmg", 3), rtols=(1e-6, 1e-4, 1e-4), + labels=labels, fingerprint=fingerprint, + what="mesh-owned FMG annulus rotated free-slip") + accuracy_anchor(values, ANCHORS["annulus_fmg"], fingerprint, labels, + what="mesh-owned FMG annulus rotated free-slip") + + def test_rotated_freeslip_spherical3d_partition_independent(): """3D spherical shell (free-slip inner+outer, all three rotation nullspace - modes): the parallel solve reproduces the serial velocity L2, converges, and + modes): the parallel solve reproduces its own np=1 velocity L2, converges, and stays within the bounded outer iteration count (the 1/mu-mass Schur preconditioner — issue #248's rotated blow-out was ~44 its).""" - L2, its, reason = _spherical3d_diagnostics() - L2_ref = GOLDEN_SPHERICAL3D + values, fingerprint, its, reason = _spherical3d_diagnostics() + # ABSOLUTE checks on the solve itself, not partition comparisons. + # + # INVARIANT: everything asserted before the `compare` below must be identical on + # every rank, because `compare` calls the COLLECTIVE `serial_reference`. A gate + # that fails on one rank only would take that rank out of the broadcast and hang + # the others instead of failing the test. `reason` and `its` come from the + # solver's own collective telemetry and `nnodes` (in the sigma test) is bcast, so + # all of them satisfy it today — but nothing enforces it, so put new absolute + # gates AFTER the compare unless you have checked. assert reason > 0, f"3D spherical rotated solve diverged: reason {reason}" assert its <= 25, f"3D spherical Schur iteration blow-out: {its} outer its" - assert np.isclose(L2, L2_ref, rtol=1e-5, atol=0), ( - f"3D spherical velocity L2 differs serial vs np={uw.mpi.size}: " - f"{L2_ref} vs {L2}") - - -@pytest.mark.xfail( - reason="#564: free-slip solves are partition dependent. CI measures a spherical " - "topography coefficient 0.4149689252074 serial vs 0.4125278837958 at np=2 " - "(5.9e-03, the largest of the family). PRE-EXISTING and not caused by " - "#560/#561: the same seven assertions fail with numbers identical to every " - "digit at #561's merge base with only the scripts/test.sh test_10*py line " - "enabled, which is how they became visible at all — this whole batch had " - "never run in CI. This case passes an explicit analytic normal=, and its " - "numbers are bit-identical before and after #560's nodal-normal fix, so it " - "is definitively not that mechanism. Passes locally on macOS/arm64, so " - "strict=False; see #564 for the full table.", - strict=False) + compare(values, _reference("spherical3d", 1), rtols=(1e-5,), + labels=("velocity L2",), fingerprint=fingerprint, + what="3D spherical rotated free-slip") + accuracy_anchor(values, ANCHORS["spherical3d"], fingerprint, ("velocity L2",), + what="3D spherical rotated free-slip") + + def test_rotated_freeslip_spherical3d_topography_partition_independent(): """3D boundary-mass recovery gives partition-independent topography coefficients.""" - coefficients = _spherical3d_topography_diagnostics() - labels = ( - "surface all", - "surface vertices", - "surface midpoints", - "CMB all", - "CMB vertices", - "CMB midpoints", - ) - for label, value, reference in zip( - labels, coefficients, GOLDEN_SPHERICAL3D_TOPO - ): - assert np.isclose(value, reference, rtol=1e-6, atol=0), ( - f"3D {label} differs serial vs np={uw.mpi.size}: " - f"{reference} vs {value}" - ) + values, fingerprint = _spherical3d_topography_diagnostics() + labels = ("surface all", "surface vertices", "surface midpoints", + "CMB all", "CMB vertices", "CMB midpoints") + compare(values, _reference("spherical3d_topo", 6), rtols=(1e-6,) * 6, + labels=labels, fingerprint=fingerprint, + what="3D spherical rotated topography") + # Zhong l=2 benchmark coefficients — physics numbers, and the reason an absolute + # anchor matters more here than anywhere else in this file. + accuracy_anchor(values, ANCHORS["spherical3d_topo"], fingerprint, labels, + what="3D spherical rotated topography") def test_rotated_freeslip_box_nonlinear_partition_independent(): """NONLINEAR rotated free-slip is partition-independent: a power-law box solved by - the manual Newton/Picard loop reproduces the serial velocity L2 and iteration count - at np=2/4 — the rotated residual/Jacobian, the increment solve and the constraint - zeroing are all parallel-safe (ownership-relative indexing, collective norms).""" - L2, iters = _box_nonlinear_diagnostics() - L2_ref, iters_ref = GOLDEN_BOX_NONLINEAR - assert np.isclose(L2, L2_ref, rtol=1e-6, atol=0), ( - f"nonlinear box velocity L2 differs serial vs np={uw.mpi.size}: {L2_ref} vs {L2}") - assert iters == iters_ref, ( - f"nonlinear iteration count differs serial vs np={uw.mpi.size}: {iters_ref} vs {iters}") + the manual Newton/Picard loop reproduces its own np=1 velocity L2 AND iteration + count at np=2/4 — the rotated residual/Jacobian, the increment solve and the + constraint zeroing are all parallel-safe (ownership-relative indexing, collective + norms).""" + values, fingerprint = _box_nonlinear_diagnostics() + # rtol=0 on the iteration count: it is an integer and the claim is exact equality. + labels = ("velocity L2", "nonlinear iteration count") + compare(values, _reference("nonlinear", 2), rtols=(1e-6, 0.0), labels=labels, + fingerprint=fingerprint, what="nonlinear box rotated free-slip") + accuracy_anchor(values, ANCHORS["nonlinear"], fingerprint, labels, + what="nonlinear box rotated free-slip") def test_rotated_freeslip_box_sigma_nn_partition_independent(): """sigma_nn (boundary_normal_traction) recovery is partition-independent: the whole- - boundary relL2 / |corr| vs analytic SolCx sigma_yy match the serial reference (and - stay accurate) in parallel — the reaction read + consistent-mass de-smear are + boundary relL2 / |corr| vs analytic SolCx sigma_yy reproduce the np=1 run (and stay + accurate) in parallel — the reaction read + consistent-mass de-smear are parallel-safe.""" - relL2, corr, nnodes = _box_sigma_diagnostics() - relL2_ref, corr_ref = GOLDEN_BOX_SIGMA + (relL2, corr, nnodes), fingerprint = _box_sigma_diagnostics() + # ABSOLUTE checks: the gathered node set, and the accuracy against the analytic + # solution. Neither is a partition comparison. assert nnodes == 49, f"expected 49 top nodes, gathered {nnodes} at np={uw.mpi.size}" - assert np.isclose(relL2, relL2_ref, rtol=1e-4, atol=0), ( - f"sigma_nn relL2 differs serial vs np={uw.mpi.size}: {relL2_ref} vs {relL2}") - assert np.isclose(corr, corr_ref, rtol=1e-4, atol=0), ( - f"sigma_nn corr differs serial vs np={uw.mpi.size}: {corr_ref} vs {corr}") assert relL2 < 0.10, f"sigma_nn relL2 vs analytic {relL2:.3f} too large" + labels = ("sigma_nn relL2", "sigma_nn |corr|") + compare((relL2, corr), _reference("sigma", 2), rtols=(1e-4, 1e-4), labels=labels, + fingerprint=fingerprint, what="box sigma_nn recovery") + accuracy_anchor((relL2, corr), ANCHORS["sigma"], fingerprint, labels, + what="box sigma_nn recovery") def test_rotated_freeslip_dynamic_topography_partition_independent(): """The dynamic_topography surface field is partition-independent: the collective - BdIntegral L2 of h over Top matches the serial reference at np=2/4. Guards the - field write (a per-node write would deadlock when a rank owns none of the boundary) - and the parallel reaction recovery underneath.""" - bdl2 = _box_topography_bdl2() - assert np.isclose(bdl2, GOLDEN_TOPO_BDL2, rtol=1e-6, atol=0), ( - f"topography BdIntegral L2 differs serial vs np={uw.mpi.size}: " - f"{GOLDEN_TOPO_BDL2} vs {bdl2}") + BdIntegral L2 of h over Top reproduces the np=1 run at np=2/4. Guards the field + write (a per-node write would deadlock when a rank owns none of the boundary) and + the parallel reaction recovery underneath.""" + values, fingerprint = _box_topography_bdl2() + compare(values, _reference("topo", 1), rtols=(1e-6,), + labels=("topography BdIntegral L2",), fingerprint=fingerprint, + what="box dynamic topography") + accuracy_anchor(values, ANCHORS["topo"], fingerprint, + ("topography BdIntegral L2",), what="box dynamic topography") + + +_DIAGNOSTICS = { + "box": _box_diagnostics, + "annulus": _annulus_diagnostics, + "annulus_fmg": _annulus_fmg_diagnostics, + "spherical3d": _spherical3d_diagnostics, + "spherical3d_topo": _spherical3d_topography_diagnostics, + "nonlinear": _box_nonlinear_diagnostics, + "sigma": _box_sigma_diagnostics, + "topo": _box_topography_bdl2, +} + + +def _reference(kind, count): + """The np=1 payload for ``kind``, with its length asserted. A reference of the + wrong length would silently shorten the ``zip`` in ``compare`` and drop + assertions.""" + reference = serial_reference(__file__, kind) + assert len(reference["values"]) >= count, ( + f"serial reference for {kind!r} has {len(reference['values'])} values, " + f"expected at least {count}") + reference = dict(reference) + reference["values"] = reference["values"][:count] + return reference if __name__ == "__main__": - # Recompute the serial GOLDEN references: - # `python {box,annulus,annulus_fmg,sigma,topo}`. + # Single-rank child of the parallel run (see serial_reference), and a + # human-readable recompute: + # `python {box,annulus,annulus_fmg,spherical3d,spherical3d_topo, + # nonlinear,sigma,topo}`. import sys _kind = sys.argv[1] if len(sys.argv) > 1 else "box" - if _kind == "nonlinear": - _L2, _its = _box_nonlinear_diagnostics() - if uw.mpi.rank == 0: - print(f"DIAG_NONLINEAR {_L2:.12e} {_its}") - elif _kind == "topo": - _b = _box_topography_bdl2() - if uw.mpi.rank == 0: - print(f"DIAG_TOPO bdl2={_b:.9e}") - elif _kind == "sigma": - _r = _box_sigma_diagnostics() - if uw.mpi.rank == 0: - print(f"DIAG_SIGMA relL2={_r[0]:.6e} corr={_r[1]:.6f} nodes={_r[2]}") - elif _kind == "annulus": - _L2, _lo, _up = _annulus_diagnostics() - if uw.mpi.rank == 0: - print(f"DIAG_ANNULUS {_L2:.12e} {_lo:.6e} {_up:.6e}") - elif _kind == "annulus_fmg": - _L2, _lo, _up = _annulus_fmg_diagnostics() - if uw.mpi.rank == 0: - print(f"DIAG_ANNULUS_FMG {_L2:.12e} {_lo:.6e} {_up:.6e}") - elif _kind == "spherical3d": - _L2, _its, _reason = _spherical3d_diagnostics() - if uw.mpi.rank == 0: - print(f"DIAG_SPHERICAL3D {_L2:.12e} its={_its} reason={_reason}") - elif _kind == "spherical3d_topo": - _cell_size = float(sys.argv[2]) if len(sys.argv) > 2 else 0.25 - _coefficients = _spherical3d_topography_diagnostics(_cell_size) - if uw.mpi.rank == 0: - print( - f"DIAG_SPHERICAL3D_TOPO cell_size={_cell_size:.8f} " - + " ".join(f"{value:.12e}" for value in _coefficients) - ) - else: - _L2, _verr = _box_diagnostics() - if uw.mpi.rank == 0: - print(f"DIAG_BOX {_L2:.12e} {_verr:.6e}") + # `python spherical3d_topo 0.3` still works: the extra argument is a + # debugging affordance (a coarser shell for a quick look) and is NOT reachable from + # serial_reference, which always calls the default so the anchor stays comparable. + _extra = [float(a) for a in sys.argv[2:]] + _result = _DIAGNOSTICS[_kind](*_extra) + # _spherical3d_diagnostics carries its solver telemetry after the fingerprint. + _values, _fingerprint = _result[0], _result[1] + emit(_values, _fingerprint) + uw.mpi.pprint(f"DIAG_{_kind.upper()} " + + " ".join(f"{v:.12e}" for v in _values) + + f" [cells={_fingerprint[0]:.0f} vol={_fingerprint[1]:.12g}]") diff --git a/tests/parallel/test_1066_rotated_datum_parallel.py b/tests/parallel/test_1066_rotated_datum_parallel.py index 652421341..ed5f6fcf0 100644 --- a/tests/parallel/test_1066_rotated_datum_parallel.py +++ b/tests/parallel/test_1066_rotated_datum_parallel.py @@ -21,26 +21,24 @@ import sympy import underworld3 as uw -pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(180)] - -# serial reference (np1), used to catch a partition-dependent datum bug -_INT_VV_REF = 6.6559607579 - - -@pytest.mark.xfail( - reason="#564: free-slip solves are partition dependent. CI measures the solve " - "energy 6.65602336 against the recorded 6.6559607579 at np=2 (9.4e-06) — " - "the smallest of the family, but the same defect. PRE-EXISTING and not " - "caused by #560/#561: the same seven assertions fail with numbers " - "identical to every digit at #561's merge base with only the " - "scripts/test.sh test_10*py line enabled, which is how they became " - "visible at all — this whole batch had never run in CI. This case passes " - "an explicit analytic normal=, and its numbers are bit-identical before " - "and after #560's nodal-normal fix, so it is definitively not that " - "mechanism. Passes locally on macOS/arm64, so strict=False; see #564 for " - "the full table.", - strict=False) -def test_rotated_datum_prescribed_normal_partition_independent(): +from serial_reference import ( + accuracy_anchor, compare, emit, mesh_fingerprint, serial_reference) + +# ABSOLUTE accuracy anchor, gated on the mesh fingerprint — `compare` proves the answer +# is partition-independent, not that it is right. The value is the pre-#568 +# `_INT_VV_REF`; the fingerprint is (owned cell count, integral 1 dV) for +# Annulus(0.5, 1.0, cellSize=0.1). rtol 1e-2, i.e. "is this still the same answer". +ANCHOR_DATUM = {"fingerprint": [600, 2.356187202481425], + "values": (6.6559607579,)} + +# The timeout covers the np=1 child the partition test spawns as well as the +# parallel solve itself. +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(900)] + + +def _datum_diagnostics(): + """Prescribe ``u.n = cos(theta)`` on the annulus surface and return + ``((relL2 of the datum residual, solve energy int v.v), fingerprint)``.""" RI, RO = 0.5, 1.0 mesh = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.1, qdegree=3) x, y = mesh.X @@ -63,12 +61,36 @@ def test_rotated_datum_prescribed_normal_partition_independent(): vn = (v.sym[0] * x + v.sym[1] * y) / r num = float(uw.maths.BdIntegral(mesh=mesh, fn=(vn - x / r) ** 2, boundary="Upper").evaluate()) den = float(uw.maths.BdIntegral(mesh=mesh, fn=(x / r) ** 2, boundary="Upper").evaluate()) - relL2 = np.sqrt(num / den) - assert relL2 < 1.0e-3, f"prescribed u.n=cos(theta) not imposed (relL2={relL2:.3e})" - + relL2 = float(np.sqrt(num / den)) vv = float(uw.maths.Integral(mesh, v.sym.dot(v.sym)).evaluate()) - assert abs(vv - _INT_VV_REF) / _INT_VV_REF < 1.0e-6, \ - f"solve energy is partition-dependent: ∫v·v={vv:.8e} vs ref {_INT_VV_REF}" + return (relL2, vv), mesh_fingerprint(mesh) + + +def test_rotated_datum_prescribed_normal_partition_independent(): + """The datum is imposed, and the solve energy reproduces THIS ENVIRONMENT's np=1 + run rather than a constant recorded elsewhere. + + The 9.4e-06 this test used to report on CI (6.65602336 measured against a recorded + 6.6559607579) was the distance to a golden taken on macOS/arm64 while the runner's + gmsh built a different annulus — not a partition effect. See the module docstring + of ``test_1064_rotated_freeslip_parallel.py`` for the CI measurement that settles + it. + """ + values, fingerprint = _datum_diagnostics() + # ABSOLUTE check: the prescribed datum really is imposed. The residual is the + # P2/faceting interpolation error, not a partition effect, so this is an accuracy + # gate and NOT a partition comparison. + assert values[0] < 1.0e-3, ( + f"prescribed u.n=cos(theta) not imposed (relL2={values[0]:.3e})") + reference = serial_reference(__file__, "datum") + compare(values[1:], {"values": reference["values"][1:], + "fingerprint": reference["fingerprint"]}, + rtols=(1e-6,), + labels=("solve energy int v.v",), fingerprint=fingerprint, + what="rotated prescribed-normal datum") + accuracy_anchor(values[1:], ANCHOR_DATUM, fingerprint, + ("solve energy int v.v",), + what="rotated prescribed-normal datum") def test_rotated_datum_nonlinear_parallel(): @@ -110,3 +132,16 @@ def test_rotated_datum_nonlinear_parallel(): den = float(uw.maths.BdIntegral(mesh=mesh, fn=(x / r) ** 2, boundary="Upper").evaluate()) relL2 = np.sqrt(max(num, 0.0) / den) assert relL2 < 1.0e-3, f"nonlinear u.n=cos(theta) not imposed (relL2={relL2:.3e})" + + +if __name__ == "__main__": + # Single-rank child of the parallel run (see serial_reference), and a + # human-readable recompute: `python datum`. + import sys + _kind = sys.argv[1] if len(sys.argv) > 1 else "datum" + if _kind != "datum": + raise SystemExit(f"unknown kind {_kind!r}") + _values, _fingerprint = _datum_diagnostics() + emit(_values, _fingerprint) + uw.mpi.pprint(f"DIAG_DATUM relL2={_values[0]:.6e} vv={_values[1]:.12e} " + f"[cells={_fingerprint[0]:.0f} vol={_fingerprint[1]:.12g}]") diff --git a/tests/parallel/test_1069_boundary_normal_parallel.py b/tests/parallel/test_1069_boundary_normal_parallel.py new file mode 100644 index 000000000..2eb43c712 --- /dev/null +++ b/tests/parallel/test_1069_boundary_normal_parallel.py @@ -0,0 +1,390 @@ +"""``mesh.boundary_normal()`` on a CURVED boundary must not depend on the partition. + +The nodal normal is ``Σ_f |f| n̂_f`` over every facet of the boundary that meets the +node. A boundary facet is labelled on exactly one rank, so a node on a partition seam +sees only SOME of its facets locally, and normalising that partial sum gives it a +rotated normal (#564). The routine used to do exactly that, under a ``TODO(parallel)`` +comment asserting it was harmless. + +It was not. On ``Annulus(cellSize=0.12)``, worst nodal normal against the exact radial +one — the pre-fix numbers these tests are calibrated to fail at: + + boundary np=1 np=2 np=3 np=4 + Upper 3.0e-10 5.8e-02 5.8e-02 5.8e-02 + Lower 5.5e-10 1.1e-01 1.1e-01 1.1e-01 + +5.8e-02 is 3.3 degrees, and it is not a rounding of the right answer: it is what you +get by taking ONE of a vertex's two facet normals instead of the average of both, so +its size is set by the facet's angular span and does not shrink with more ranks. +``mesh.boundary_normal()`` is the default constraint direction for +``add_constraint_bc`` and ``add_nitsche_bc``, and that error moved a constrained +free-slip answer by 3.4 % between np=1 and np=2 (#495, one member of #564). + +Flat axis-aligned walls are immune — every facet of such a wall carries the same +normal, so a partial stencil normalises to the same answer — which is why the box +tests never saw this and why every test here is on a curved boundary or a corner. + +AND IT WAS ALSO WRONG IN SERIAL, IN 3-D +--------------------------------------- +Facet contributions used to reach their DOFs by a kd-tree query for "the nodes +nearest the facet centroid". On a TETRAHEDRAL boundary the three DOFs nearest a face +centroid are not always that face's own three vertices, so the query picked up a +neighbour and the assembled normal was wrong on a **uniform** mesh **at np=1**: + + SphericalShell(0.55, 1.0, cs=0.35) vs the global facet sum + Upper old 4.71e-02 new 1.9e-16 + Lower old 1.03e-01 new 2.2e-16 (1.03e-01 is 5.9 degrees) + +2-D is unaffected (annulus 1.1e-16 old and new; box bit-identical) — on an edge the +two nearest DOFs to the midpoint are always its own two vertices. + +That defect is invisible to the analytic radial oracle: at this resolution the honest +faceting error is 8.8e-02 / 2.26e-01 and the WORST node is the same before and after, +so ``max‖n − r̂‖`` is identical either way. Only the global-facet-sum oracle sees it. +That is why this file has one, and why **this file deliberately carries no +``mpi(min_size=2)`` mark**: it must run at np=1, or the serial 3-D path — which no +other test in the tree reaches, every serial default-normal test being on a box — +goes uncovered again. The single test that genuinely needs np>1 skips itself. + +Run with: + python -m pytest tests/parallel/test_1069_boundary_normal_parallel.py + mpirun -n 2 python -m pytest --with-mpi \\ + tests/parallel/test_1069_boundary_normal_parallel.py + mpirun -n 4 python -m pytest --with-mpi \\ + tests/parallel/test_1069_boundary_normal_parallel.py +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from underworld3.discretisation.discretisation_mesh import Mesh +from underworld3.utilities.facet_normals import facet_measure_and_normal + +from serial_reference import compare, emit, mesh_fingerprint, serial_reference + +# NO mpi(min_size=2): see "AND IT WAS ALSO WRONG IN SERIAL, IN 3-D" above. The one +# test that cannot mean anything at np=1 skips itself. +pytestmark = [pytest.mark.timeout(600)] + +_KEY_DECIMALS = 10 + + +def _key(coord): + return tuple(np.round(np.asarray(coord, dtype=float).ravel(), _KEY_DECIMALS)) + + +def _assembled_normals(mesh, boundary): + """``{coordinate key: normal}`` for every node the assembled field gave a normal, + merged across ranks, plus the worst disagreement between two ranks holding the + same node. + + The disagreement is expected to be exactly zero even on the UNFIXED code: writing + the field through ``var.data`` runs a local→global (INSERT) → global→local round + trip, which overwrites every ghost copy with the owner's value. Ranks agreeing is + therefore NOT evidence that the normal is right — the owner can be confidently + wrong and everyone will copy it. The load-bearing assertions below are the two + oracles, not this number; it is checked because a non-zero value would mean the + sync assumption has changed underneath us. + """ + mesh.boundary_normal(boundary) + var = mesh._boundary_normal_vars[boundary] + coords = np.asarray(var.coords) + data = np.asarray(var.data) + live = np.linalg.norm(data, axis=1) > 0.5 + local = {_key(c): tuple(float(t) for t in v) + for c, v in zip(coords[live], data[live])} + + merged, disagreement = {}, 0.0 + for rank_values in uw.mpi.comm.allgather(local): + for key, value in rank_values.items(): + if key in merged: + disagreement = max(disagreement, max( + abs(a - b) for a, b in zip(merged[key], value))) + else: + merged[key] = value + return merged, disagreement + + +def _global_facet_sum_oracle(mesh, boundary): + """``{coordinate key: normal}`` computed the way the routine SHOULD compute it, + by an independent route: gather every boundary facet in the mesh onto every rank, + de-duplicate by centroid, and do the whole ``Σ_f |f| n̂_f`` in numpy. + + This is an oracle rather than a re-run of the code under test: it never touches the + variable's section, its local↔global scatter, or the DM at all after the gather, so + it cannot reproduce a plumbing bug in the reduction. It is partition-independent by + construction — the facet set it sums over is the global one on every rank. + """ + dm = mesh.dm + cdim = mesh.cdim + vertex_start, vertex_end = dm.getDepthStratum(0) + coord_section = dm.getCoordinateSection() + coord_vector = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, cdim) + + local_facets = {} + for facet in mesh._boundary_facets(boundary): + measure, normal, exterior = facet_measure_and_normal(dm, facet) + if not exterior: + continue + _, centroid, _ = dm.computeCellGeometryFVM(facet) + vertices = [ + _key(coord_vector[coord_section.getOffset(point) // cdim]) + for point in dm.getTransitiveClosure(facet)[0] + if vertex_start <= int(point) < vertex_end + ] + local_facets[_key(centroid)] = (measure, tuple(normal[:cdim]), vertices) + + all_facets = {} + for rank_facets in uw.mpi.comm.allgather(local_facets): + all_facets.update(rank_facets) + + accumulated = {} + for measure, normal, vertices in all_facets.values(): + for vertex in vertices: + accumulated[vertex] = accumulated.get(vertex, np.zeros(cdim)) \ + + measure * np.asarray(normal) + return {k: tuple(v / np.linalg.norm(v)) for k, v in accumulated.items()} + + +def _worst_difference(left, right): + """Largest ``|a - b|`` over the keys the two dictionaries share, and the number of + keys only one of them has (which must be zero — a node set that depends on the + partition is its own defect).""" + shared = set(left) & set(right) + assert shared, "the two normal fields share no nodes — the comparison is vacuous" + worst = max(float(np.linalg.norm(np.asarray(left[k]) - np.asarray(right[k]))) + for k in shared) + return worst, len(set(left) ^ set(right)) + + +def _annulus(cell_size=0.12): + return uw.meshing.Annulus(radiusInner=0.5, radiusOuter=1.0, + cellSize=cell_size, qdegree=3) + + +# --------------------------------------------------------------------------- # +# Oracle 1 — the exact radial normal on an annulus +# --------------------------------------------------------------------------- # + +@pytest.mark.parametrize("boundary,sign", [("Upper", 1.0), ("Lower", -1.0)]) +def test_boundary_normal_annulus_matches_exact_radial(boundary, sign): + """On a circular arc the measure-weighted average of a vertex's two chord normals + is EXACTLY radial by symmetry, so the exact normal is an oracle with no golden and + no serial run: serial lands at 3.0e-10 (Upper) / 5.5e-10 (Lower), and any partial + stencil lands at ~half the facet's angular span, 5.8e-02 / 1.1e-01. The 1e-6 gate + sits four orders below the defect and four above the discretisation floor. + + ``sign`` is the DOMAIN's outward direction: +r̂ on the outer arc, −r̂ on the inner + one (the #560 convention — outward is away from the facet's own support cell, which + on a concave boundary points at the centre of curvature). + """ + mesh = _annulus() + normals, disagreement = _assembled_normals(mesh, boundary) + assert normals, f"no nodes carry a normal on {boundary}" + worst = max( + float(np.linalg.norm( + np.asarray(n) - sign * np.asarray(k) / np.linalg.norm(k))) + for k, n in normals.items()) + assert worst < 1.0e-6, ( + f"{boundary} nodal normal is {worst:.3e} from the exact radial one at " + f"np={uw.mpi.size} over {len(normals)} nodes — a partial (rank-local) facet " + f"stencil, #564") + assert disagreement == 0.0, ( + f"{boundary} normals differ by {disagreement:.3e} between ranks holding the " + f"same node") + + +# --------------------------------------------------------------------------- # +# Oracle 2 — the global facet sum, in 2-D and 3-D +# --------------------------------------------------------------------------- # + +def _shell(cell_size=0.35): + return uw.meshing.SphericalShell(radiusInner=0.55, radiusOuter=1.0, + cellSize=cell_size, qdegree=3) + + +def test_boundary_normal_shell_matches_global_facet_sum_in_serial(): + """THE SERIAL 3-D CASE, on its own, because it is a second defect this change + fixes and nothing else in the tree reaches it. + + At np=1 there is no partition, so the cross-rank reduction is a no-op and the ONLY + thing under test is where a facet's contribution lands. The old kd-tree route + ("the DOFs nearest the facet centroid") put it on the wrong vertex often enough to + be 4.71e-02 (Upper) and 1.03e-01 (Lower, ≈5.9°) away from the true facet sum on a + UNIFORM shell; the section-based route is 1.9e-16 / 2.2e-16. + + Kept separate from the parametrised test below, and asserted at np=1 explicitly, + so that this cannot quietly become parallel-only again — which is exactly how the + defect survived: `test_1069` was `mpi(min_size=2)`, and every other serial test of + the default normal is on a box, where flat walls make the question vacuous. + """ + if uw.mpi.size != 1: + pytest.skip("this assertion is about the SERIAL path; the parametrised " + "facet-sum test covers np>1") + mesh = _shell() + for boundary in ("Upper", "Lower"): + assembled, _ = _assembled_normals(mesh, boundary) + worst, missing = _worst_difference( + assembled, _global_facet_sum_oracle(mesh, boundary)) + assert missing == 0, f"shell {boundary}: {missing} nodes in only one node set" + assert worst < 1.0e-12, ( + f"shell {boundary} at np=1: assembled normal is {worst:.3e} from the " + f"global facet sum over {len(assembled)} nodes — the facet-to-DOF routing " + f"is wrong, independently of any partition effect") + + +@pytest.mark.parametrize("geometry", ["annulus", "shell"]) +def test_boundary_normal_matches_global_facet_sum(geometry): + """The assembled field equals the sum over the GLOBAL facet set, in 2-D and 3-D. + + This is the assertion that generalises: a 3-D shell's faceted normal is 8.8e-02 + (outer) / 2.3e-01 (inner) from the exact radial one at this resolution — that is + honest discretisation error, not a defect, so the analytic oracle above cannot be + used there. What must hold is that every rank count produces the sum over ALL the + facets, which is what this compares against. + """ + mesh = _annulus() if geometry == "annulus" else _shell() + for boundary in ("Upper", "Lower"): + assembled, _ = _assembled_normals(mesh, boundary) + oracle = _global_facet_sum_oracle(mesh, boundary) + worst, missing = _worst_difference(assembled, oracle) + assert missing == 0, ( + f"{geometry} {boundary}: {missing} nodes are in one of the assembled / " + f"oracle node sets and not the other at np={uw.mpi.size}") + assert worst < 1.0e-12, ( + f"{geometry} {boundary}: assembled normal differs from the global facet " + f"sum by {worst:.3e} at np={uw.mpi.size} over {len(assembled)} nodes " + f"(#564: a rank-local stencil)") + + +def test_boundary_normal_oracle_fires_without_the_reduction(): + """NEGATIVE CONTROL. With the cross-rank completion disabled — which is exactly the + pre-#564 code — both oracles above must FAIL. Without this, a fix that quietly + stopped the assembly from producing anything, or an oracle that re-derived the same + wrong answer, would pass the suite. + """ + if uw.mpi.size == 1: + pytest.skip("a rank-local stencil is the complete stencil in serial") + + original = Mesh._sum_local_dofs_across_ranks + Mesh._sum_local_dofs_across_ranks = lambda self, subdm, values: values + try: + mesh = _annulus() + normals, _ = _assembled_normals(mesh, "Upper") + radial_error = max( + float(np.linalg.norm( + np.asarray(n) - np.asarray(k) / np.linalg.norm(k))) + for k, n in normals.items()) + oracle_error, _ = _worst_difference( + normals, _global_facet_sum_oracle(mesh, "Upper")) + finally: + Mesh._sum_local_dofs_across_ranks = original + + assert radial_error > 1.0e-3, ( + f"the radial oracle does not see a rank-local stencil at np={uw.mpi.size} " + f"(error {radial_error:.3e}); it cannot be trusted to see a regression") + assert oracle_error > 1.0e-3, ( + f"the global-facet-sum oracle does not see a rank-local stencil at " + f"np={uw.mpi.size} (error {oracle_error:.3e})") + + +# --------------------------------------------------------------------------- # +# The regression the fix is most likely to introduce +# --------------------------------------------------------------------------- # + +def test_boundary_normal_corner_is_not_averaged_across_boundaries(): + """Each boundary is assembled into its OWN variable, so the vertex where Top meets + Right keeps (0,1) for Top and (1,0) for Right rather than the 45-degree bisector. + A cross-rank reduction done on a shared vector instead of per boundary would + average across the discontinuity and this is what would catch it.""" + mesh = uw.meshing.StructuredQuadBox( + elementRes=(8, 8), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=3) + top, _ = _assembled_normals(mesh, "Top") + right, _ = _assembled_normals(mesh, "Right") + corner = _key((1.0, 1.0)) + assert corner in top and corner in right, ( + "the (1,1) corner should carry a normal on BOTH Top and Right") + assert np.allclose(top[corner], (0.0, 1.0), rtol=0, atol=1e-12), ( + f"Top normal at the corner is {top[corner]}, not (0, 1)") + assert np.allclose(right[corner], (1.0, 0.0), rtol=0, atol=1e-12), ( + f"Right normal at the corner is {right[corner]}, not (1, 0)") + + +# --------------------------------------------------------------------------- # +# End to end: the two consumers of the default normal +# --------------------------------------------------------------------------- # + +def _nitsche_annulus_diagnostics(): + """Nitsche free-slip on a curved boundary through the DEFAULT (assembled) normal — + ``add_nitsche_bc`` shares ``mesh.boundary_normal`` with ``add_constraint_bc`` and + was equally exposed, with no test covering it. Returns ((velocity L2, radial + leakage on Upper), mesh fingerprint). + + With the cross-rank completion disabled this reads, at solver tolerance 1e-9: + + np=1 L2 1.413526513414e-02 leak 1.729537946e-04 + np=2 1.414677409720e-02 3.752196269e-04 + np=3 1.422846013917e-02 2.873963117e-04 + np=4 1.413111891962e-02 3.149424657e-04 + + — 6.6e-03 in the velocity and more than DOUBLE the wall-normal leakage. With it, + every rank count agrees to 3.6e-10 in the velocity and bit-identically in the + leakage. Both are stable from tolerance 1e-9 to 1e-12, so neither is the linear + solve. + + ``local_h=False`` is deliberate and it is not a workaround for this fix. The + default ``local_h=True`` scales the Nitsche penalty by ``mesh.cell_size()``, which + is built from ``Mesh._get_mesh_sizes`` — a kd-tree query against THIS RANK's cell + centroids, and so partition-dependent in its own right (on this mesh the field's + sum is 26.0822 at np=1, 26.1211 at np=2, 26.1386 at np=4, and its max moves at + np=4). That is a SEPARATE defect from the boundary normal, it is not what #564 is + about, and leaving it in would make this test measure the two together. See the + TODO(BUG) on ``Mesh._assemble_cell_size``. + """ + RI, RO = 0.5, 1.0 + mesh = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.12, qdegree=3) + x, y = mesh.X + r = sympy.sqrt(x ** 2 + y ** 2) + theta = sympy.atan2(y, x) + v = uw.discretisation.MeshVariable("Vn69", mesh, mesh.dim, degree=2, continuous=True) + p = uw.discretisation.MeshVariable("Pn69", mesh, 1, degree=1, continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.bodyforce = sympy.Matrix([[ + x / r * sympy.cos(4 * theta) * (r - RI) * (RO - r) * 40.0, + y / r * sympy.cos(4 * theta) * (r - RI) * (RO - r) * 40.0]]) + stokes.add_essential_bc((0.0, 0.0), "Lower") + # default normal= is the assembled one — that is what is under test + stokes.add_nitsche_bc(0.0, "Upper", local_h=False) + stokes.tolerance = 1.0e-9 + stokes.petsc_options["snes_type"] = "ksponly" + stokes.solve() + + L2 = float(np.sqrt(uw.maths.Integral(mesh, v.sym.dot(v.sym)).evaluate())) + vr = v.sym[0] * x / r + v.sym[1] * y / r + leak = float(np.sqrt(uw.maths.BdIntegral( + mesh=mesh, fn=vr ** 2, boundary="Upper").evaluate())) + return (L2, leak), mesh_fingerprint(mesh) + + +def test_nitsche_freeslip_annulus_partition_independent(): + """``add_nitsche_bc`` on a curved boundary with the default normal reproduces its + OWN np=1 answer (computed in this environment, not a recorded constant).""" + values, fingerprint = _nitsche_annulus_diagnostics() + compare(values, serial_reference(__file__, "nitsche"), + rtols=(1e-8, 1e-8), labels=("velocity L2", "Upper radial leakage"), + fingerprint=fingerprint, what="Nitsche free-slip annulus") + + +if __name__ == "__main__": + import sys + + _kind = sys.argv[1] if len(sys.argv) > 1 else "nitsche" + if _kind == "nitsche": + _values, _fingerprint = _nitsche_annulus_diagnostics() + emit(_values, _fingerprint) + uw.mpi.pprint(f"DIAG_NITSCHE L2={_values[0]:.12e} leak={_values[1]:.6e}") + else: + raise SystemExit(f"unknown kind {_kind!r}")