From 4d660447ea45810f6af61e2543f3c6c7e3663616 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 11:10:38 +1000 Subject: [PATCH 1/9] A boundary normal built from one rank's facets is a rotated normal: complete the sum across ranks `mesh.boundary_normal()` is the default constraint direction for both `add_constraint_bc` and `add_nitsche_bc`. It accumulated `|f| * n_f` over THIS RANK's labelled boundary facets and then normalised. 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 a partial sum does not give a slightly noisy normal - it gives the wrong one. Measured on Annulus(cellSize=0.12), the worst nodal normal against the exact radial one: 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. It does not shrink with more ranks because it is not round-off: 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. End to end that moved a constrained free-slip velocity by 3.4% between np=1 and np=2 - #495's number to every digit, and one of the seven rows of #564. The source carried a TODO(parallel) prescribing exactly this reduction and asserting the code was "parallel-SAFE as-is (only the handful of partition-seam surface vertices get a slightly-rotated unit normal)". That claim was the bug; it is replaced by the measurement. The weighted contributions are now summed through the variable's OWN sub-DM local-to-global scatter before normalising, so every rank ends with a bit-identical normal at a shared node. The reduction is COLLECTIVE - a rank owning no facet of this boundary still takes part. A plain ADD is exact with no de-duplication because no boundary facet is labelled twice and none is labelled away from its owner: measured at np=2,3,4 on the annulus and the spherical shell, for the per-boundary label this uses AND for the consolidated UW_Boundaries label. Work vectors are created from the sub-DM rather than borrowed from the variable, which is what the TODO's abandoned first attempt SEGV'd on. Two things went with it. The DOF a facet contributes to now comes from the variable's own section on that sub-DM instead of a kd-tree lookup of the DOFs nearest the facet centroid: the kd-tree is a heuristic that can mis-assign on a graded mesh, and it cannot be made to agree across ranks because each rank's tree is built from its own local coordinates - which would have made the sum land on different points on different ranks. And the outward-orientation rule, the facet measure and the internal-boundary guard now come from ONE shared `utilities/facet_normals.facet_measure_and_normal`, used by all three copies of this accumulation (rotated_bc, this one, and boundary_flux._node_normals, whose copy was still on the pre-#560 rule that orients against the mean of the rank's coordinates and points INTO the domain on a concave boundary). Three copies drifting is what produced #564 twelve weeks after #560 fixed the first one. Serial results are bit-identical: 3.005063e-10 / 5.488346e-10 before and after, and now the same to every digit at np=1, 2, 3 and 4. The new tests carry their own oracles rather than a recorded constant: on a circular arc the measure-weighted average of a vertex's two chord normals is EXACTLY radial by symmetry, and in 3-D - where faceting is honest discretisation error, 8.8e-02 on a coarse shell - the assembled field is compared against the sum over the GLOBAL facet set, gathered and done in numpy. The negative control disables the reduction and confirms both oracles fire at np>1. Corner preservation is asserted on a box, because averaging across a boundary discontinuity is the regression this change could most easily introduce. While measuring the Nitsche twin, a SECOND and unrelated partition dependence turned up: `mesh.cell_size()`, the default Nitsche local h, is built from a kd-tree query against this rank's centroids and moves with the rank count (field sum 26.0822 / 26.1211 / 26.1386 at np=1/2/4). It is marked TODO(BUG) where it lives and the Nitsche test passes local_h=False so it measures the normal and not the two together. Not fixed here - `_get_mesh_sizes` also feeds get_min_radius, the adaptivity metrics and the free-surface relaxation. Fixes the constrained half of #564 (and #495). The rotated rows of #564 are a separate question. Underworld development team with AI support from Claude Code --- docs/developer/subsystems/rotated-freeslip.md | 39 ++ .../discretisation/discretisation_mesh.py | 164 ++++++--- src/underworld3/utilities/boundary_flux.py | 29 +- src/underworld3/utilities/facet_normals.py | 64 ++++ src/underworld3/utilities/rotated_bc.py | 60 ++-- tests/parallel/serial_reference.py | 133 +++++++ .../test_1069_boundary_normal_parallel.py | 335 ++++++++++++++++++ 7 files changed, 734 insertions(+), 90 deletions(-) create mode 100644 src/underworld3/utilities/facet_normals.py create mode 100644 tests/parallel/serial_reference.py create mode 100644 tests/parallel/test_1069_boundary_normal_parallel.py 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/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 03fe527f2..780dce6df 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,12 @@ 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. + Returns a sympy Matrix (row) of the P1 normal-field components, for use as the constraint direction in Nitsche/penalty BCs. @@ -2975,56 +2981,115 @@ 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``. + """ + 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)) + ncomp = var.num_components + # 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)) face_pts = self._boundary_facets(name) - 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. + # 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 — a heuristic that can mis-assign on a graded mesh, and one + # that 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: + ssec = subdm.getLocalSection() + for f in face_pts: + # 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] + + if dm.comm.getSize() > 1: + accum = self._sum_local_dofs_across_ranks(subdm, accum) + finally: + indexset.destroy() + subdm.destroy() + 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. + """ + 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() + # A DOF constrained out of the global vector would come back zero; keep + # this rank's own contribution there rather than losing it. The mesh DM + # carries no essential-BC constraints (those live on the solver DM), so + # this does not fire today — it is here so that it degrades to the old + # rank-local answer instead of to zero if that ever changes. + dropped = (~summed.any(axis=1)) & values.any(axis=1) + summed[dropped] = values[dropped] + return summed + def cell_size(self): """Local, per-cell characteristic mesh size as a scalar field symbol. @@ -3106,6 +3171,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: 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..dfdb8fe76 --- /dev/null +++ b/src/underworld3/utilities/facet_normals.py @@ -0,0 +1,64 @@ +"""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): + """``(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, in which case the normal has been oriented AWAY from + that cell — the domain's outward direction. When it is False the normal carries + PETSc's own sign and no orientation claim is made. + + 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). + """ + 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: + _, 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/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..ddc6f0e18 --- /dev/null +++ b/tests/parallel/serial_reference.py @@ -0,0 +1,133 @@ +"""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=1800): + """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.""" + 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): + env = {k: v for k, v in os.environ.items() + if not k.startswith(_MPI_ENV_PREFIXES)} + try: + 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 {os.path.basename(module_file)}:{kind} " + f"timed out after {timeout}s") + for line in proc.stdout.splitlines(): + if line.startswith("SERIALREF "): + return json.loads(line[len("SERIALREF "):]) + return (f"serial reference for {os.path.basename(module_file)}:{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 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.""" + 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)}") + fp_note = (f" [mesh np=1: cells={ref_fp[0]:.0f} vol={ref_fp[1]:.12g}; " + f"np={uw.mpi.size}: cells={fingerprint[0]:.0f} " + f"vol={fingerprint[1]:.12g}]") + 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}") 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..af2178801 --- /dev/null +++ b/tests/parallel/test_1069_boundary_normal_parallel.py @@ -0,0 +1,335 @@ +"""``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. + +Run with: + 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 + +pytestmark = [pytest.mark.mpi(min_size=2), 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) + + +@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}") From 91384dd3a59e2f11d27dc389b3f63fae01261c88 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 11:25:44 +1000 Subject: [PATCH 2/9] The constrained partition tests now agree at every rank count: drop their xfails test_1063's three assertions were xfail(strict=False) against #564. With the boundary normal completed across ranks they hold on their own terms, so the xfails go. Measured here, this file's own diagnostics: [iso] velocity L2 np=1 6.194547793939e-01 np=2 6.194547793947e-01 np=3 6.194547793476e-01 np=4 6.194547793462e-01 [ti] velocity L2 np=1 3.925981603048e-01 np=2 3.925981604290e-01 np=3 3.925981603382e-01 np=4 3.925981604662e-01 - spreads of 1.2e-10 and 4.1e-10. Before the fix the same runs gave [iso] 5.982807168537e-01 at np=2 (-3.4%, and #495's number to every digit) and 6.187878061677e-01 at np=4, with the topography 2.4% adrift; [ti] moved 0.34% at np=2 and 2.8% at np=4. Reproduced by disabling only the cross-rank sum, so the attribution is a bisection and not an inference. The reference is no longer a constant recorded on a developer's machine. It is this file's own np=1 run, computed in the same environment by a single-rank child (tests/parallel/serial_reference.py), which is the property the test claims to be testing. A stored constant additionally measures the host's mesh generator: four other rows of #564 turned out to be exactly that, and the message "differs serial vs np=N" on a comparison against a number from somewhere else is what sent that investigation the wrong way for a day. Both mesh fingerprints (owned cell count, integral 1 dV) are printed on failure so a mesh difference reads as a mesh difference. Two tolerances are stated rather than inherited. Velocity keeps 1e-8, the parallel reduction order. Topography moves to 1e-5 from 1e-6 because it is read off the multiplier, whose [p,lambda] Schur sub-block grinds into its 200-iteration cap on this problem - the ti case sits at 6.4e-07 against the old 1e-6 gate, which is a 1.6x margin and not a tolerance, while 1e-5 still leaves three orders below the 2.4% it has to catch. The raw mean pressure is now asserted ABSOLUTELY (|meanP| < 1e-6) instead of against the reference: it is pinned to zero by the automatic gauge, and comparing two numbers near machine zero relatively measures nothing. Negative controls, both run: with the cross-rank sum reverted, [iso] and [ti] fail at np=2; with a 1e-3 relative partition dependence injected into every diagnostic, all three fail. Neither passes vacuously. Underworld development team with AI support from Claude Code --- ...test_1063_constrained_freeslip_parallel.py | 168 +++++++++--------- 1 file changed, 80 insertions(+), 88 deletions(-) diff --git a/tests/parallel/test_1063_constrained_freeslip_parallel.py b/tests/parallel/test_1063_constrained_freeslip_parallel.py index 843f99a08..0060154b1 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,22 @@ import underworld3 as uw -pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(180)] +from serial_reference import compare, emit, mesh_fingerprint, serial_reference + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(600)] -# 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), -} +# 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 +80,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 +92,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 +135,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 +145,42 @@ 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}") + compare((L2, topo), {"values": [reference["values"][0], reference["values"][2]], + "fingerprint": reference["fingerprint"]}, + rtols=(_RTOL_VELOCITY, _RTOL_TOPOGRAPHY), + labels=("velocity L2", "mean-stripped topography"), + fingerprint=fingerprint, 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) + compare(values, serial_reference(__file__, kind), + rtols=(_RTOL_VELOCITY, _RTOL_TOPOGRAPHY), + labels=("velocity L2", "mean-stripped topography"), + fingerprint=fingerprint, 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)) From 5575f8862584dd894080b0edec0f198f527f260d Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 11:26:04 +1000 Subject: [PATCH 3/9] The rotated partition tests were comparing against another host: make them self-referential The four remaining #564 xfails are not a defect in the solver. Running test_1064's own _annulus_diagnostics on ONE CI host at both rank counts gives np=1 1.897329151623790e-02 and np=2 1.897329151623740e-02 - agreement to the 13th significant figure, with leak_lo and leak_up identical to every digit - while BOTH differ from the recorded GOLDEN_ANNULUS by exactly the same +1.676e-04. The rotated free-slip path is partition-independent. gmsh builds a different triangulation on the Linux runner than on macOS/arm64, the goldens were recorded on macOS, and because these files are mpi(min_size=2) CI has never run np=1 to notice. So the defect is in the test design, and it is the kind that costs real time: the assertion said "differs serial vs np=N" while what it measured was "differs from a number recorded elsewhere". Those are different statements. Seven assertions failing under the first reading looked like one family of solver bugs; four of them were the mesh. Every partition assertion in test_1064 and test_1066 now compares against a np=1 run of the same file, computed in the same environment by a single-rank child (tests/parallel/serial_reference.py, which scrubs the launcher's MPI variables so the child runs as a singleton rather than trying to join the parent job). Whatever mesh the host generates, both sides of the comparison use it. Each failure message carries both runs' mesh fingerprints - owned cell count and integral 1 dV, both partition-independent - so a future host difference is legible instead of mysterious. The checks that are genuinely ABSOLUTE are kept and relabelled as accuracy checks, with a comment saying they are NOT partition-independence checks: the analytic SolCx velocity error, sigma_nn's relL2 against the exact solution, the gathered top-node count, the Schur iteration ceiling, the converged reason, and the datum residual. Conflating the two is what this commit exists to stop. The mesh-owned FMG test deliberately takes the EXPLICIT-registration run as its np=1 reference: #467 is precisely the claim that the two routes produce the same solve, so that comparison asserts both properties at once. All four xfails removed. 11 passed at np=2 and at np=4. Negative control: with a 1e-3 relative partition dependence injected into every diagnostic, all eleven partition assertions in the two files fire. This addresses the ROTATED rows of #564 as a test-design problem. It does not touch the solver, and it is a separate matter from the constrained rows, which are fixed by the boundary-normal reduction earlier in this branch. Underworld development team with AI support from Claude Code --- .../test_1064_rotated_freeslip_parallel.py | 361 ++++++++---------- .../test_1066_rotated_datum_parallel.py | 74 ++-- 2 files changed, 200 insertions(+), 235 deletions(-) diff --git a/tests/parallel/test_1064_rotated_freeslip_parallel.py b/tests/parallel/test_1064_rotated_freeslip_parallel.py index 63a935a11..baa148ed2 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,11 @@ 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 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)] def _wrap(dm, m0): @@ -105,7 +95,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 +127,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 +157,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 +224,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 +276,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 +307,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 +353,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 +378,48 @@ 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. + assert values[1] < 1e-3, ( + f"box velocity error {values[1]:.2e} too large at np={uw.mpi.size}") + + 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}") + compare(values, _reference("annulus", 3), rtols=(1e-6, 1e-4, 1e-4), + labels=("velocity L2", "Lower radial leakage", "Upper radial leakage"), + fingerprint=fingerprint, 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() + compare(values, _reference("annulus_fmg", 3), rtols=(1e-6, 1e-4, 1e-4), + labels=("velocity L2", "Lower radial leakage", "Upper radial leakage"), + fingerprint=fingerprint, what="custom-FMG annulus rotated free-slip") def test_rotated_freeslip_mesh_owned_fmg_pickup(): @@ -456,151 +431,117 @@ 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) + compare(values, _reference("annulus_fmg", 3), rtols=(1e-6, 1e-4, 1e-4), + labels=("velocity L2", "Lower radial leakage", "Upper radial leakage"), + fingerprint=fingerprint, 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. 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") + + 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") 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. + compare(values, _reference("nonlinear", 2), rtols=(1e-6, 0.0), + labels=("velocity L2", "nonlinear iteration count"), + fingerprint=fingerprint, 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" + compare((relL2, corr), _reference("sigma", 2), rtols=(1e-4, 1e-4), + labels=("sigma_nn relL2", "sigma_nn |corr|"), + fingerprint=fingerprint, 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") + + +_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}") + _result = _DIAGNOSTICS[_kind]() + # _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..0234f0336 100644 --- a/tests/parallel/test_1066_rotated_datum_parallel.py +++ b/tests/parallel/test_1066_rotated_datum_parallel.py @@ -21,26 +21,16 @@ 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 compare, emit, mesh_fingerprint, serial_reference + +# 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 +53,33 @@ 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") def test_rotated_datum_nonlinear_parallel(): @@ -110,3 +121,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}]") From a4dc58c0eb59e4231b1a5bba4a1c7122eab13653 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 11:58:07 +1000 Subject: [PATCH 4/9] Keep the new collective out of reach of an existing swallow-and-continue Completing the boundary-normal sum across ranks turned _assemble_boundary_normal into a COLLECTIVE routine, and deform() calls it inside a `try: ... except: pass` that predates the change. A rank-local failure there - the case that comment names, a boundary whose label vanished from the current DM after region extraction - would now take one rank out of the reduction and HANG the ranks that did find the label, converting a swallowed error into a deadlock. The facet walk is the only rank-local part, so its guard moves INSIDE the collective: a rank that cannot complete it contributes zero and still takes part in the reduction. The outer guard in deform() can then only ever see a symmetric failure, and says so. Also corrects scripts/test.sh's comment to name test_1069 and #564 alongside #560 in the parallel solver batch it enables. Underworld development team with AI support from Claude Code --- scripts/test.sh | 7 +-- .../discretisation/discretisation_mesh.py | 51 ++++++++++++------- 2 files changed, 38 insertions(+), 20 deletions(-) diff --git a/scripts/test.sh b/scripts/test.sh index 569ab4558..afa0dc223 100755 --- a/scripts/test.sh +++ b/scripts/test.sh @@ -118,9 +118,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 780dce6df..05690d5d1 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -3016,8 +3016,6 @@ def _assemble_boundary_normal(self, var, name): # it is simply a row that stays zero until the reduction fills it. accum = numpy.zeros_like(numpy.asarray(var.data)) - face_pts = self._boundary_facets(name) - # 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 @@ -3027,22 +3025,33 @@ def _assemble_boundary_normal(self, var, name): indexset, subdm = dm.createSubDM(var.field_id) try: ssec = subdm.getLocalSection() - for f in face_pts: - # 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: + try: + 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 - accum[ssec.getOffset(q) // ncomp] += measure * nrm[:cdim] + 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: + # Sanctioned swallow, and it is deliberately INSIDE the collective: + # the facet walk is rank-local, and a rank that cannot complete it + # (a boundary whose label vanished from the current DM, e.g. after + # region extraction) simply contributes nothing. Raising or + # returning here would take this rank out of the reduction below + # and HANG the ranks that did find the label — the failure mode the + # caller's `except: pass` in deform() would otherwise convert a + # one-rank error into. + accum[...] = 0.0 if dm.comm.getSize() > 1: accum = self._sum_local_dofs_across_ranks(subdm, accum) @@ -3987,6 +3996,14 @@ def _do_move(): # 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. + # + # NB _assemble_boundary_normal is COLLECTIVE (it completes + # the facet sum across ranks, #564). It swallows its own + # rank-local failures INSIDE that collective so one rank + # cannot drop out of the reduction and hang the others; this + # outer guard therefore only ever sees a symmetric failure. + # The dict is built by a collective accessor, so every rank + # iterates the same boundaries in the same order. pass # Likewise refresh the local cell-size field (Nitsche penalty scaling) # so its cell-constant data tracks the deformed geometry. From bb3c2819a460f063eef4438fc58373bd21428fc8 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 13:35:22 +1000 Subject: [PATCH 5/9] The fault normal was a fourth copy of the same six lines, not a different rule The PR excluded fault_contact._fault_pair_nodes from the shared helper on the grounds that it "walks an internal surface where both support cells exist by design and the orientation IS the plus/minus side split - a different rule, not a stale copy". Review checked that against the source and it does not hold. The loop is the SAME six lines as the block extracted from rotated_bc: same computeCellGeometryFVM, same unit-normalise with the identical +1e-30, same flip against getSupport(f)[0]'s centroid, same float(vol) measure weight, same full-length dot. The one difference is that the getSupportSize == 1 guard is dropped so the flip is unconditional. That is a different GUARD, not a different rule - and it leaves a fourth verbatim copy of the six lines whose duplication this PR itself argues is the root cause pattern. So it consumes the helper. `orient="support0"` names the fault's rule explicitly: flip away from support[0] unconditionally, which on a split fault is the plus/minus side split and is coherent along the surface because every facet of the surface lists the same side first. The default `orient="exterior"` is the boundary rule - orient only where "outward from the domain" is defined. The docstring says which surfaces each is safe on and why the wrong choice makes neighbouring facets CANCEL in a measure-weighted sum. Behaviour-preserving by construction (the two branches are line-for-line what the callers had) and by measurement: tests/test_0846_fault_contact.py, test_0845_fault_split.py, test_0847_fault_api.py, test_0850_faults.py - 48 passed. Also records in the helper docstring that the orientation dot product is taken over dm.getCoordinateDim() components, which matters only for a manifold mesh where a caller slices the result to fewer. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/facet_normals.py | 31 ++++++++++++++++++---- src/underworld3/utilities/fault_contact.py | 22 ++++++++++----- 2 files changed, 41 insertions(+), 12 deletions(-) diff --git a/src/underworld3/utilities/facet_normals.py b/src/underworld3/utilities/facet_normals.py index dfdb8fe76..057543e6b 100644 --- a/src/underworld3/utilities/facet_normals.py +++ b/src/underworld3/utilities/facet_normals.py @@ -39,25 +39,46 @@ import numpy as np -def facet_measure_and_normal(dm, facet): +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, in which case the normal has been oriented AWAY from - that cell — the domain's outward direction. When it is False the normal carries - PETSc's own sign and no orientation claim is made. + 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: + 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 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 From 9270aca4dc40c032981ab9a197a6273af85d169c Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 13:35:22 +1000 Subject: [PATCH 6/9] A rank-local failure must be loud and collective, not a silent zero normal or a hang Review of #568 found the new error path was worse than the defect it guards. Both halves measured, both fixed, both with a negative control. SILENT ZERO NORMALS. The `except Exception: accum[...] = 0.0` caught anything - including an IndexError from the new row arithmetic - and turned it into "this rank contributed nothing". Measured at np=2 with rank 1's facet walk raising: deform() returned cleanly on both ranks, nothing was printed, and rank 1's OWNED boundary DOFs came back with a ZERO normal (sum|n| 2.186 against a control of 21.306). The constraint direction over that part of the boundary is then the zero vector, with a converged solve. That is strictly worse than the 3.3 degrees this branch exists to remove, and indistinguishable from success. HANG ON AN ASYMMETRIC RAISE. The previous guard covered only the facet walk. Everything before `dm.createSubDM` - var.num_components, the zeros_like on var.data - and `subdm.getLocalSection()` sat outside it, and createSubDM is collective. Measured: rank 1 walked out through deform()'s `except: pass` and rank 0 blocked in the sub-DM collective until the launcher killed the job. One discipline applied twice. Every rank-local step now sets a failure flag instead of escaping; the flag is ALL-REDUCED before the reduction so every rank takes the same branch; then every rank raises, carrying the failing ranks' messages. deform() decides collectively too - its all-reduce sits on the exception path, so "every rank takes the same branch" is a property of that loop rather than an assumption inherited from the callee - and it WARNS instead of skipping in silence, because a swallowed refresh leaves a stale constraint direction on a moved boundary. Negative controls at np=2, both green: a failure inside the facet walk, and a failure before createSubDM, each raise on BOTH ranks with the originating rank's message, both reach the following barrier, and neither leaves a zero normal. A third injection that replaces the ENTIRE routine still hangs and cannot be defended against from the caller - rank 0 is inside createSubDM before any flag could be exchanged. What remains unguarded inside the routine is only createSubDM itself, the two destroys, the allgather on the raise path and the final collective pack: all collective, so all symmetric. A rank-local raise added anywhere else would reintroduce the hang, which is why the docstring says so. Drops the "the round trip lost this DOF, keep the local value" fallback in _sum_local_dofs_across_ranks. It asked "did it come back all-zero?" when the question is "was it constrained out?", and those differ on exactly the input that matters: a node whose GLOBAL contributions cancel would have its rank-local PARTIAL value restored - #564 again, on the one case the reduction exists to get right. Dead code, therefore also untested and unreachable by the negative control. getConstraintDof is the predicate if it is ever needed. Refuses a non-degree-1 variable rather than half-filling it: the row arithmetic assumes one node per DMPlex point, and boundary_normal() adopts a pre-existing `_n_bd_` variable if one is registered, so a checkpoint restore or user code could hand it a P3 space whose edge carries two nodes. And discloses the SERIAL 3-D correction this branch makes - see the following commit for the coverage. Underworld development team with AI support from Claude Code --- .../discretisation/discretisation_mesh.py | 177 ++++++++++++++---- 1 file changed, 138 insertions(+), 39 deletions(-) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 05690d5d1..ecaf3eb1c 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2923,6 +2923,25 @@ def boundary_normal(self, boundary): (#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. @@ -3005,27 +3024,74 @@ def _assemble_boundary_normal(self, var, name): (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 - ncomp = var.num_components - # 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)) + 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 — a heuristic that can mis-assign on a graded mesh, and one - # that cannot be made to agree across ranks because each rank's tree is - # built from its own local coordinates.) + # 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: - ssec = subdm.getLocalSection() 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 @@ -3042,23 +3108,29 @@ def _assemble_boundary_normal(self, var, name): if ssec.getDof(q) <= 0: continue accum[ssec.getOffset(q) // ncomp] += measure * nrm[:cdim] - except Exception: - # Sanctioned swallow, and it is deliberately INSIDE the collective: - # the facet walk is rank-local, and a rank that cannot complete it - # (a boundary whose label vanished from the current DM, e.g. after - # region extraction) simply contributes nothing. Raising or - # returning here would take this rank out of the reduction below - # and HANG the ranks that did find the label — the failure mode the - # caller's `except: pass` in deform() would otherwise convert a - # one-rank error into. - accum[...] = 0.0 - - if dm.comm.getSize() > 1: - accum = self._sum_local_dofs_across_ranks(subdm, accum) + 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() + 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])}") + mag = numpy.sqrt(numpy.sum(accum ** 2, axis=1)) nonzero = mag > 1.0e-30 accum[nonzero] /= mag[nonzero, numpy.newaxis] @@ -3077,6 +3149,18 @@ def _sum_local_dofs_across_ranks(self, subdm, values): 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() @@ -3090,13 +3174,6 @@ def _sum_local_dofs_across_ranks(self, subdm, values): finally: lvec.destroy() gvec.destroy() - # A DOF constrained out of the global vector would come back zero; keep - # this rank's own contribution there rather than losing it. The mesh DM - # carries no essential-BC constraints (those live on the solver DM), so - # this does not fire today — it is here so that it degrades to the old - # rank-local answer instead of to zero if that ever changes. - dropped = (~summed.any(axis=1)) & values.any(axis=1) - summed[dropped] = values[dropped] return summed def cell_size(self): @@ -3986,25 +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. + # 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. # - # NB _assemble_boundary_normal is COLLECTIVE (it completes - # the facet sum across ranks, #564). It swallows its own - # rank-local failures INSIDE that collective so one rank - # cannot drop out of the reduction and hang the others; this - # outer guard therefore only ever sees a symmetric failure. - # The dict is built by a collective accessor, so every rank - # iterates the same boundaries in the same order. - pass + # 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: From a2c7a5cdfb30a74f374bf5f2a40ba14261b25d7f Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 13:35:59 +1000 Subject: [PATCH 7/9] This also fixes a SERIAL 3-D defect: say so, and cover it at np=1 Undisclosed in the first version of this branch, and the most valuable thing in it. 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, where there is no partition effect at all: SphericalShell(0.55, 1.0, cs=0.35), np=1, vs the global facet sum Upper old 4.712169e-02 new 1.922963e-16 Lower old 1.029976e-01 new 2.225049e-16 (1.03e-01 is 5.9 degrees) 2-D is untouched: annulus 1.241267e-16 old and new on Upper, box bit-identical - on an edge the two nearest DOFs to the midpoint are always its own two vertices. Reproduced here by running the pre-branch assembly verbatim alongside the new one in a single serial process, so the only variable is the DOF routing. The branch said "serial results are bit-identical". That was true of what had been measured (2-D) and false in 3-D, so mesh.boundary_normal() - the default constraint direction of add_constraint_bc and add_nitsche_bc, and of FreeSurface's wall-normal datum - has been returning materially wrong normals on every 3-D curved boundary in serial, and this branch corrects them. That belongs in the docstring and in the release notes, not in a diff. NOTHING COVERED IT, and the shape of the miss is worth recording: the analytic radial oracle CANNOT see it. 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 - rhat| is identical either way. Only the global-facet-sum oracle catches it - and that oracle lived in a file marked mpi(min_size=2), so it never ran at np=1. Every other serial test of the default normal is on a box, where flat walls make the question vacuous. The mark is gone. test_1069 now runs at every rank count - 7 passed, 1 skipped in 8.7 s at np=1, with the one genuinely MPI-only test (the reduction's negative control) skipping itself through a guard the file already carried; 21 passed, 1 skipped at np=2, 3 and 4. A dedicated serial 3-D assertion is added rather than relying on the parametrised one, so this cannot quietly become parallel-only again, and scripts/test.sh runs the file in the SERIAL job as well - none of the serial batches reached this path before, and test_106*py, where most of the default-normal consumers live, is deliberately not batched at all. Underworld development team with AI support from Claude Code --- scripts/test.sh | 8 +++ .../test_1069_boundary_normal_parallel.py | 57 ++++++++++++++++++- 2 files changed, 64 insertions(+), 1 deletion(-) diff --git a/scripts/test.sh b/scripts/test.sh index afa0dc223..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. diff --git a/tests/parallel/test_1069_boundary_normal_parallel.py b/tests/parallel/test_1069_boundary_normal_parallel.py index af2178801..2eb43c712 100644 --- a/tests/parallel/test_1069_boundary_normal_parallel.py +++ b/tests/parallel/test_1069_boundary_normal_parallel.py @@ -24,7 +24,30 @@ 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 \\ @@ -40,7 +63,9 @@ from serial_reference import compare, emit, mesh_fingerprint, serial_reference -pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(600)] +# 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 @@ -179,6 +204,36 @@ def _shell(cell_size=0.35): 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. From 79c348b2f9ff013bbbd43b79f6fb47aa7b98f5bc Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 13:35:59 +1000 Subject: [PATCH 8/9] Restore the eleven absolute anchors behind a mesh-fingerprint gate Making the partition tests self-referential was right and it removed something that was doing real work. Self-reference 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 - each passes a self-referential test. Eleven constants were deleted and only some were replaced by an accuracy check, so nothing pinned the annulus radial leakages, the custom-FMG answer, the Zhong coefficients, test_1066's energy, or test_1063's velocity and topography. The reason for deleting them was real - they are host-specific, because gmsh triangulates differently across platforms, and a mismatch then reads as a physics regression. But that is an argument for gating them, which is what the #564 investigation actually recommended, and this branch had already built mesh_fingerprint() and then not used it that way. So every golden comes back as an explicitly-labelled ACCURACY anchor at rtol=1e-2, in front of which accuracy_anchor() checks the mesh: same triangulation, the gate is live; different, it SKIPS with both fingerprints in the message instead of failing for the wrong reason. Partition independence is asserted either way, three to eight orders tighter. All eleven values are the pre-branch goldens, verified to reproduce on this host: iso 6.194547793955e-01, ti 3.925981604039e-01, annulus leakages 4.563841e-05 / 9.341699e-06, FMG 1.906961759626e-02, spherical3d 4.069689334228e-03, the six Zhong coefficients, nonlinear (8.069396188270e-04, 8), sigma (5.554578e-02, 0.998466), topo 2.553916470e-01, datum 6.6559607579. compare() now ASSERTS the two fingerprints match rather than only printing them. Its docstring promised "the day that stops being true it says so", which was only kept if the test failed for another reason first. Two more from the same review: * serial_reference deadlocked on any child-launch failure except a timeout. An OSError from subprocess.run, a truncated SERIALREF line, a MemoryError on capture_output all propagated out of rank 0 BEFORE the broadcast, leaving every other rank in MPI_Bcast - busy-polling a core each - until pytest's own timeout fired. _run_child now catches everything and returns the failure as the string the contract already had for timeouts. Its default timeout drops 1800 -> 600 s so it is shorter than the mark that wraps it and can actually fire. * the absolute gates that run BEFORE the collective serial_reference() call (nnodes, its, reason) are all collective or bcast today, so they are symmetric - but nothing enforced it and a one-rank failure there would hang rather than fail. Stated as an invariant where the next person will read it. Also restores the cell_size debugging argument to test_1064's __main__, which the dict dispatch had dropped. Underworld development team with AI support from Claude Code --- tests/parallel/serial_reference.py | 108 +++++++++++++++--- ...test_1063_constrained_freeslip_parallel.py | 40 ++++++- .../test_1064_rotated_freeslip_parallel.py | 102 +++++++++++++++-- .../test_1066_rotated_datum_parallel.py | 13 ++- 4 files changed, 231 insertions(+), 32 deletions(-) diff --git a/tests/parallel/serial_reference.py b/tests/parallel/serial_reference.py index ddc6f0e18..ae147c3a2 100644 --- a/tests/parallel/serial_reference.py +++ b/tests/parallel/serial_reference.py @@ -71,10 +71,23 @@ def mesh_fingerprint(mesh): return [cells, volume] -def serial_reference(module_file, kind, timeout=1800): +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.""" + 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] @@ -88,19 +101,28 @@ def serial_reference(module_file, kind, timeout=1800): def _run_child(module_file, kind, timeout): - env = {k: v for k, v in os.environ.items() - if not k.startswith(_MPI_ENV_PREFIXES)} + """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 {os.path.basename(module_file)}:{kind} " - f"timed out after {timeout}s") - for line in proc.stdout.splitlines(): - if line.startswith("SERIALREF "): - return json.loads(line[len("SERIALREF "):]) - return (f"serial reference for {os.path.basename(module_file)}:{kind} printed no " + 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:]}") @@ -115,19 +137,77 @@ def emit(values, fingerprint): "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.""" + 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)}") - fp_note = (f" [mesh np=1: cells={ref_fp[0]:.0f} vol={ref_fp[1]:.12g}; " - f"np={uw.mpi.size}: cells={fingerprint[0]:.0f} " - f"vol={fingerprint[1]:.12g}]") + 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 0060154b1..260ee4684 100644 --- a/tests/parallel/test_1063_constrained_freeslip_parallel.py +++ b/tests/parallel/test_1063_constrained_freeslip_parallel.py @@ -45,10 +45,36 @@ import underworld3 as uw -from serial_reference import compare, emit, mesh_fingerprint, serial_reference +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 @@ -152,11 +178,13 @@ def test_constrained_raw_gauge_partition_independent(): # 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=("velocity L2", "mean-stripped topography"), + 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"]) @@ -165,10 +193,12 @@ def test_constrained_freeslip_partition_independent(kind): 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=("velocity L2", "mean-stripped topography"), + 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__": diff --git a/tests/parallel/test_1064_rotated_freeslip_parallel.py b/tests/parallel/test_1064_rotated_freeslip_parallel.py index baa148ed2..46f8f05d4 100644 --- a/tests/parallel/test_1064_rotated_freeslip_parallel.py +++ b/tests/parallel/test_1064_rotated_freeslip_parallel.py @@ -57,12 +57,54 @@ from underworld3.function import analytic as A from underworld3.utilities import custom_mg -from serial_reference import compare, emit, mesh_fingerprint, serial_reference +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): """Wrap a refined DMPlex as a UW3 mesh carrying the base mesh's boundaries and @@ -389,6 +431,8 @@ def test_rotated_freeslip_box_partition_independent(): values, fingerprint = _box_diagnostics() compare(values[:1], _reference("box", 1), rtols=(1e-8,), labels=("velocity L2",), fingerprint=fingerprint, what="box rotated free-slip") + accuracy_anchor(values[:1], ANCHORS["box"], fingerprint, ("velocity L2",), + 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. assert values[1] < 1e-3, ( @@ -407,9 +451,13 @@ def test_rotated_freeslip_annulus_partition_independent(): 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. + labels = ("velocity L2", "Lower radial leakage", "Upper radial leakage") compare(values, _reference("annulus", 3), rtols=(1e-6, 1e-4, 1e-4), - labels=("velocity L2", "Lower radial leakage", "Upper radial leakage"), - fingerprint=fingerprint, what="annulus rotated free-slip") + 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(): @@ -417,9 +465,12 @@ def test_rotated_freeslip_annulus_fmg_partition_independent(): 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.""" 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=("velocity L2", "Lower radial leakage", "Upper radial leakage"), - fingerprint=fingerprint, what="custom-FMG annulus rotated free-slip") + 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(): @@ -438,9 +489,12 @@ def test_rotated_freeslip_mesh_owned_fmg_pickup(): 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=("velocity L2", "Lower radial leakage", "Upper radial leakage"), - fingerprint=fingerprint, what="mesh-owned FMG annulus rotated free-slip") + 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(): @@ -450,11 +504,21 @@ def test_rotated_freeslip_spherical3d_partition_independent(): preconditioner — issue #248's rotated blow-out was ~44 its).""" 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" 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(): @@ -465,6 +529,10 @@ def test_rotated_freeslip_spherical3d_topography_partition_independent(): 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(): @@ -475,9 +543,11 @@ def test_rotated_freeslip_box_nonlinear_partition_independent(): norms).""" values, fingerprint = _box_nonlinear_diagnostics() # rtol=0 on the iteration count: it is an integer and the claim is exact equality. - compare(values, _reference("nonlinear", 2), rtols=(1e-6, 0.0), - labels=("velocity L2", "nonlinear iteration count"), + 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(): @@ -490,9 +560,11 @@ def test_rotated_freeslip_box_sigma_nn_partition_independent(): # solution. Neither is a partition comparison. assert nnodes == 49, f"expected 49 top nodes, gathered {nnodes} at np={uw.mpi.size}" assert relL2 < 0.10, f"sigma_nn relL2 vs analytic {relL2:.3f} too large" - compare((relL2, corr), _reference("sigma", 2), rtols=(1e-4, 1e-4), - labels=("sigma_nn relL2", "sigma_nn |corr|"), + 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(): @@ -504,6 +576,8 @@ def test_rotated_freeslip_dynamic_topography_partition_independent(): 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 = { @@ -538,7 +612,11 @@ def _reference(kind, count): # nonlinear,sigma,topo}`. import sys _kind = sys.argv[1] if len(sys.argv) > 1 else "box" - _result = _DIAGNOSTICS[_kind]() + # `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) diff --git a/tests/parallel/test_1066_rotated_datum_parallel.py b/tests/parallel/test_1066_rotated_datum_parallel.py index 0234f0336..ed5f6fcf0 100644 --- a/tests/parallel/test_1066_rotated_datum_parallel.py +++ b/tests/parallel/test_1066_rotated_datum_parallel.py @@ -21,7 +21,15 @@ import sympy import underworld3 as uw -from serial_reference import compare, emit, mesh_fingerprint, serial_reference +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. @@ -80,6 +88,9 @@ def test_rotated_datum_prescribed_normal_partition_independent(): 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(): From eee97812a22d71c44ea665902abb10537f6cf569 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Sat, 15 Aug 2026 14:50:28 +1000 Subject: [PATCH 9/9] Assert the accuracy gate before the anchor that can skip it The SolCx velocity-error assertion sat after accuracy_anchor(), which SKIPS on a host whose mesh fingerprint differs - so on such a host the accuracy check was dropped along with the anchor. Moved above it, with the reason stated. Also corrects the PR body: TWO tolerances moved in test_1063, not one - the velocity gate 1e-9 -> 1e-8 as well as the topography gate, and the velocity move is forced by this PR's own measurement (a [ti] spread of 4.1e-10 on 0.3926 is 1.05e-9 relative, which fails the old gate). Underworld development team with AI support from Claude Code --- tests/parallel/test_1064_rotated_freeslip_parallel.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/tests/parallel/test_1064_rotated_freeslip_parallel.py b/tests/parallel/test_1064_rotated_freeslip_parallel.py index 46f8f05d4..aeed50d4d 100644 --- a/tests/parallel/test_1064_rotated_freeslip_parallel.py +++ b/tests/parallel/test_1064_rotated_freeslip_parallel.py @@ -431,12 +431,14 @@ def test_rotated_freeslip_box_partition_independent(): values, fingerprint = _box_diagnostics() compare(values[:1], _reference("box", 1), rtols=(1e-8,), labels=("velocity L2",), fingerprint=fingerprint, what="box rotated free-slip") - accuracy_anchor(values[:1], ANCHORS["box"], fingerprint, ("velocity L2",), - 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():