diff --git a/.claude/skills/adapt-on-top-faults/SKILL.md b/.claude/skills/adapt-on-top-faults/SKILL.md index 48cd277f..db7ff2cc 100644 --- a/.claude/skills/adapt-on-top-faults/SKILL.md +++ b/.claude/skills/adapt-on-top-faults/SKILL.md @@ -190,6 +190,10 @@ isotropic viscosity). Rotated strong free-slip imposes `u·n̂=0` as an ESSENTIA constraint in a per-node (n,t) frame → machine-zero leakage AND composes with TI. ```python +# normal=None (the default) is measure-weighted and consistent with the assembly — +# prefer it. An analytic nhat is exact for the TRUE circle but keeps a consistency +# error against the faceted integral (#560); use it only when the constraint must +# follow the geometry rather than the mesh. nhat = mesh.CoordinateSystem.unit_e_0 # exact radial normal (annulus/sphere) stokes.add_rotated_freeslip_bc(0, "Upper", normal=nhat) stokes.add_rotated_freeslip_bc(0, "Lower", normal=nhat) diff --git a/CLAUDE.md b/CLAUDE.md index b5f471dd..0f7868a1 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -330,8 +330,12 @@ The PETSc-based solvers are carefully optimized and validated. **NO CHANGES with to impose `v·n̂ = 0`: - Enforces zero wall-normal flow to **machine precision** (Nitsche / penalty leak ~1e-3). -- Correct on **curved / tilted / deformed** boundaries — the normal is taken per node - (pass an analytic `normal`, e.g. `X/|X|`, for an exact normal on curved faces). +- Correct on **curved / tilted / deformed** boundaries — the normal is taken per node, + measure-weighted so it matches the straight-facet integral the assembler evaluates + (#560). Leave `normal=None` unless the constraint must follow the TRUE surface rather + than the mesh; an analytic `normal` (e.g. `X/|X|`) is exact for the geometry but keeps + a consistency error against the faceted assembly. See + `docs/developer/subsystems/rotated-freeslip.md` ("Which normal to use"). - Works **inside the nonlinear SNES** and with **geometric FMG**. It honours `solver.consistent_jacobian`: use `True` (consistent Newton) for smooth nonlinear rheologies; `"continuation"` (staged Picard→Newton) for robustness far from the diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md index 7c47271e..f1ed4cc1 100644 --- a/docs/developer/subsystems/rotated-freeslip.md +++ b/docs/developer/subsystems/rotated-freeslip.md @@ -16,14 +16,98 @@ stokes.add_rotated_freeslip_bc(0, "Upper", normal=nhat) # free-slip stokes.add_rotated_freeslip_bc(h_dot.sym[0], "Upper", normal=nhat) # u·n̂ = field ``` -`normal=None` uses the geometric facet normal; a sympy `1×dim` matrix in -`mesh.X` supplies an analytic normal (exact `X/|X|` on curved boundaries — the -preferred choice there); a constant array is also accepted. The datum must be a +`normal=None` uses the geometric facet normal (the default, and the one +consistent with what the assembler integrates — see below); a sympy `1×dim` +matrix in `mesh.X` supplies an analytic normal, exact for the TRUE surface +(`X/|X|` on a spherical cap, a constant on a planar face); a constant array is +also accepted. The datum must be a *scalar* (a number, an expression of `mesh.X`, or a scalar field read); on an enclosed boundary it must be discretely flux-free for incompressibility. A corner or 3D-edge node shared between rotated boundaries has no single normal and stays at the free-slip pinning (the datum is ignored there). +### The nodal normal is measure-weighted, not a bisector + +A node that sits on more than one facet — a vertex in 2D, a vertex or an +edge-midpoint in 3D — gets one nodal normal, while the assembler integrates the +boundary term facet by facet. The two only agree when the node's normal is +parallel to the **measure-weighted** sum `Σ_f |f| n̂_f` (edge length in 2D, face +area in 3D), which is what the geometric path accumulates. Plain bisector +averaging `Σ_f n̂_f` — what UW3 did before issue #560 — is right only where the +facets are equal; on a **kinked** wall with unequal facets it leaves a residual +`sin(Δ/2)·(|f₁|−|f₂|)/6` in the node's free tangential row (Δ = kink angle), the +exact constant-pressure vector stops being a null vector of the constrained +operator, and the pressure gauge goes unpinned. **Axis-aligned** walls are +unchanged to the last bit — their facet normals have exactly 0/±1 components, so +`Σ_f |f| n̂_f` normalises to the same floats as `Σ_f n̂_f` whatever the weights. +A flat but *tilted* wall is not covered by that argument and can move by one ulp. + +The sum runs over ALL facets meeting the node, so it must be completed **across +ranks**. Each boundary facet is labelled on exactly one rank, so a node on a +partition seam sees only some of its facets locally; the contributions are +summed through the DM's local↔global scatter before normalising, which is what +makes the normal partition-independent. Two things had to go with it: the +outward test now points away from the facet's own support cell (the mean of the +rank's coordinates is rank-local, and would let two facets of one node cancel), +and the node list comes from the local mesh's exterior facets rather than the +labelled subset, because a rank can own a node whose labelled facets are all on +neighbours. + +### Which way is "outward" — and the sign of σ_nn on an inner boundary + +The geometric normal points away from the facet's own support cell, which is the +**domain's** outward normal on any boundary. On a concave boundary — an annulus +or spherical-shell **inner** arc, the CMB — that points *toward* the centre of +curvature. + +This changed at #560. The old rule pointed away from the mean of the mesh +coordinates, which on an inner arc is *into* the domain. So on a concave +boundary, through the geometric normal: + +| quantity | before #560 | after | +|---|---:|---:| +| nodal radial component, annulus `Lower` | +1.000000 | **−1.000000** | +| `boundary_normal_traction("Lower")` | −5.233110e-02 | **+5.233110e-02** | + +Magnitudes are identical to every digit; only the sign moves. `dynamic_topography_field` +is `h = −σ_nn/(Δρ g)` on top of that number, so it reverses there too, as does the +sign of a non-zero prescribed wall-normal datum (`u·n̂ = ũ_n`: positive now means +outflow *from the domain* on an inner arc, where before it meant inflow). Convex +boundaries — every box wall, an outer arc, a spherical cap — are unaffected: the two +rules agree there. + +An **analytic** `normal=` is applied exactly as supplied and is *not* reoriented. +So `X/|X|` on an inner arc is inward-of-domain and gives σ_nn of the opposite sign +to the default. That is deliberate — the override means "use exactly this +direction", and silently flipping it would change the meaning of a user's datum — +but it means **you must pass `-X/|X|` on an inner boundary if you want the +domain-outward convention.** `test_1018_rotated_nodal_normal.py` pins both halves. + +### Which normal to use + +They answer different questions, and the trade is measurable. `|A z|/|A|_F` on +an annulus (`cellSize=0.15`), `z` = the attached constant pressure: + +| boundary | np | geometric (default) | analytic `X/\|X\|` | +|---|---:|---:|---:| +| uniform arcs | 1 | **6.6e-20** | 1.1e-14 | +| uniform arcs | 4 | **6.7e-20** | 1.1e-14 | +| skewed (non-uniform facets) | 1 | **6.6e-20** | 3.1e-07 | +| skewed (non-uniform facets) | 4 | **6.7e-20** | 3.1e-07 | + +The geometric normal is *consistent with the assembly*: it is the direction the +straight-facet boundary integral actually sees, so the constant pressure stays a +null vector to machine precision at every rank count. The analytic normal is +*consistent with the geometry*: it is tangent to the true surface, which the +faceted mesh only approximates — so the assembler and the constraint disagree by +an amount that grows with facet non-uniformity, and #560 does not remove it (the +analytic column is unchanged by this fix, and identical at every rank count). + +Prefer the default. Reach for `normal=` when the constraint must follow the true +surface rather than the mesh — a coarse spherical shell where faceting, not the +gauge, is the dominant error — and be aware that the pressure gauge is then only +as good as the numbers above. + Why strong rather than Nitsche/penalty: the constraint holds to machine precision (a penalty leaks ~1e-3, and the leak grows exactly where anisotropy makes the boundary condition matter), it is correct on curved/tilted/deformed diff --git a/scripts/test.sh b/scripts/test.sh index 37ccfb3a..569ab455 100755 --- a/scripts/test.sh +++ b/scripts/test.sh @@ -117,9 +117,12 @@ if [ $PARALLEL_RANKS -gt 0 ]; then echo "Testing global statistics and parallel operations..." mpirun -n $PARALLEL_RANKS python -m pytest --with-mpi tests/parallel/test_075*py || status=1 - # Add other parallel test categories as they're created: - # echo "Testing parallel solvers..." - # mpirun -n $PARALLEL_RANKS python -m pytest --with-mpi tests/parallel/test_10*py || status=1 + # Parallel SOLVER tests. This line was commented out, so test_1017 and + # test_1062..test_1068 — the whole rotated / constrained / MG parallel set, + # including the partition-independence guard for the rotated nodal normal + # (#560) — executed at NO rank count in CI. + echo "Testing parallel solvers..." + mpirun -n $PARALLEL_RANKS python -m pytest --with-mpi tests/parallel/test_10*py || status=1 # echo "Testing parallel I/O..." # mpirun -n $PARALLEL_RANKS python -m pytest --with-mpi tests/parallel/test_io*py || status=1 diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index d20ad05f..3ff11af0 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -6138,11 +6138,16 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): boundary : str Boundary label to constrain. normal : None or sympy 1×dim Matrix or array, optional - Per-node outward normal source. ``None`` uses the geometric facet - normal (PETSc ``computeCellGeometryFVM``; works in 2D and 3D). A - sympy ``1×dim`` matrix supplies an analytic normal (exact - ``X/|X|`` on a spherical cap, a constant on a planar face) — preferred - on curved boundaries. A constant array is also accepted. + Per-node outward normal source. ``None`` (the default, and normally + the right choice) uses the geometric facet normal, measure-weighted + so that it is consistent with the straight-facet boundary integral + the assembler evaluates. A sympy ``1×dim`` matrix supplies an + analytic normal (``X/|X|`` on a spherical cap, a constant on a + planar face): exact for the TRUE surface, but the assembler still + integrates over the facets, so it keeps a consistency error that + grows with facet non-uniformity. Use it when the constraint must + follow the geometry rather than the mesh. A constant array is also + accepted. See ``docs/developer/subsystems/rotated-freeslip.md``. Notes ----- diff --git a/src/underworld3/systems/free_surface.py b/src/underworld3/systems/free_surface.py index 75a7c030..6fd13177 100644 --- a/src/underworld3/systems/free_surface.py +++ b/src/underworld3/systems/free_surface.py @@ -279,6 +279,14 @@ def _normal_direction(self, coords): """Per-node unit vectors along the topography direction that the surface increment is deformed along — vertical (Cartesian) or radial (annulus / spherical shell); dimension-general.""" + # TODO(BUG): this is unconditionally +radial, but since #560 the geometric + # constraint normal (and so the recovered σ_nn and h_∞) is the DOMAIN's + # outward normal, which on a CONCAVE surface — an inner arc / CMB free + # surface — is −radial. The relaxation would then drive the surface away + # from equilibrium instead of toward it. Unreachable today: every curved + # free surface in the repo passes an explicit normal=rhat on an OUTER + # boundary, and every Cartesian one is a flat Top. Needs the sign taken from + # the same source as h_∞ before an inner free surface is supported. if self._radial: r = np.linalg.norm(coords, axis=1) r[r == 0.0] = 1.0 diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index 0f3cfd10..aad5153f 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/src/underworld3/utilities/boundary_flux.py @@ -213,6 +213,13 @@ def _node_normals(solver, boundary, normal, nodes, dm, dim, cvec, csec, v0, v1): """Per-node outward unit normal (only needed to project a vector reaction). ``normal`` is None (geometric facet normal), a sympy 1×dim Matrix (analytic, lambdified), or a constant (dim,) vector.""" + # TODO(BUG): the geometric branch below is a stale copy of the pre-#560 rule. + # It orients against the mean of the mesh coordinates, which is rank-local (it + # averages only this rank's points) and points INTO the domain on a concave + # boundary — rotated_bc._boundary_velocity_nodes now orients away from the + # facet's own support cell and sums across ranks. Currently unreachable: the + # only caller guards it with `if normal is not None`, so the geometric branch + # never runs. It will be wrong the day someone wires a geometric normal in here. interior_ref = cvec.mean(axis=0) sym_fn = const = None if normal is not None: diff --git a/src/underworld3/utilities/fault_contact.py b/src/underworld3/utilities/fault_contact.py index 4b0c13f3..1d3e54f8 100644 --- a/src/underworld3/utilities/fault_contact.py +++ b/src/underworld3/utilities/fault_contact.py @@ -548,9 +548,19 @@ def _fault_pair_nodes(solver, boundary): dm.getLabel(plus_name).getStratumIS(value).getIndices() if fS <= int(p) < fE] + # Facet normals accumulated to the pair nodes, weighted by the facet MEASURE for + # the same reason as the wall normals in rotated_bc (#560): the assembler + # integrates facet by facet, so a node on two facets is only consistent when its + # normal is parallel to Σ_f |f| n̂_f. The fault is not protected by being an + # interior surface — a constant pressure cancels exactly in the MEAN rows (the two + # sides carry opposite outward normals) but DOUBLES in the jump rows, and the + # jump-tangential (slip) row is free, so the bisector leaves √2·p·sin(δ/2)·(|f₁|−|f₂|)/6 + # there: a pressure-driven spurious slip at every kink node, plus the same lost + # pressure gauge. Rank-local by construction — a seam-touching fault is + # redistributed onto one rank before the split, so no cross-rank sum is needed. nacc = {} for f in facets: - _, cent, nrm = dm.computeCellGeometryFVM(f) + vol, cent, nrm = dm.computeCellGeometryFVM(f) ne = np.asarray(nrm, dtype=float) ne = ne / (np.linalg.norm(ne) + 1e-30) support = dm.getSupport(f) @@ -559,7 +569,7 @@ def _fault_pair_nodes(solver, boundary): ne = -ne for q in (int(c) for c in dm.getTransitiveClosure(f)[0]): if lsec.getFieldDof(q, _VELOCITY_FIELD) > 0: - nacc[q] = nacc.get(q, np.zeros(dim)) + ne + nacc[q] = nacc.get(q, np.zeros(dim)) + float(vol) * ne override = _compiled_normal_override(solver, boundary) if override is not None: diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 0882f854..47dfed76 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -104,14 +104,71 @@ def _boundary_velocity_nodes(solver, boundary, normal=None): unit normal. Dimension-general (2D edges, 3D faces). `normal` selects the normal source (per boundary): - * None — geometric facet normal from PETSc ``computeCellGeometryFVM`` - (area-weighted, accumulated to the facet closure points). + * None — geometric facet normal from PETSc ``computeCellGeometryFVM``, + accumulated to the facet closure points weighted by the facet + MEASURE (edge length in 2D, face area in 3D). The weight is + what makes the constraint consistent with the assembly — see + the note below. * sympy 1×dim Matrix — an analytic normal (function of mesh.X); evaluated - at each node's coordinate. Best for exact curved/planar faces - (radial ``X/|X|`` on a spherical cap; a constant on a planar - side of a regional spherical box). + at each node's coordinate (radial ``X/|X|`` on a spherical + cap; a constant on a planar side of a regional spherical + box). Exact for the TRUE surface, which is a different thing + from being consistent with the straight-facet assembly — see + "Which normal to use" in + ``docs/developer/subsystems/rotated-freeslip.md``. * (dim,) array — a constant normal vector. Returns ``[(point, n̂), ...]``. + + Why the geometric normals are MEASURE-weighted (issue #560) + ----------------------------------------------------------- + A boundary node sitting on more than one facet — a vertex in 2D, a vertex or an + edge-midpoint in 3D — gets ONE nodal normal, while the assembler integrates the + boundary term facet by facet. The rotated frame's free tangential row therefore + collects, for a constant pressure :math:`p`, + + .. math:: + r_i = p \\sum_f \\left(\\int_f \\phi_i\\, ds\\right)\\, (\\hat t_i \\cdot n_f) + + which vanishes for every tangent :math:`\\hat t_i` only when the nodal normal is + parallel to :math:`\\sum_f (\\int_f \\phi_i\\, ds)\\, n_f`. For simplicial P1 and P2 + velocity the basis integral over a facet is the facet measure times a constant that + does not depend on which facet it is, so the measure IS the right weight. That is + EXACT for every simplicial facet (2D P2 vertex :math:`|f|/6`; 3D P1 vertex + :math:`|f|/3`; 3D P2 edge-midpoint :math:`|f|/3`; the 3D P2 vertex integral is + identically zero, so that row is consistent whatever normal it is given) and for + every 2D facet at any degree, a facet being a 1D element there. On a **non-affine + 3D quad facet** (a deformed hex) the Jacobian varies across the facet, so the + measure is the leading-order weight rather than the exact one — measured, a + deformed hex box retains a residual ~7-16x its flat control while the bisector + leaves 3e8-2e9 times more, so the weighting is a large improvement there but not + exact. + + Normalising each facet normal to unit length before accumulating — what this + function used to do — gives the BISECTOR :math:`n_1 + n_2` instead, which is + correct only where the facets are equal. On a kinked boundary with unequal facets + the constant-pressure vector then stops being a null vector of the constrained + operator, the pressure gauge goes unpinned, and the answer picks up a round-off + seeded offset (#560). + + An axis-aligned wall is unchanged to the last bit: its facet normals have exactly + 0/±1 components, so :math:`\\sum_f |f| n_f` normalises to the same floats as + :math:`\\sum_f n_f` whatever the weights. A flat but TILTED wall can move by one + ulp, which is why this reads "axis-aligned" and not "flat". + + Parallel + -------- + The sum is over ALL facets meeting the node, so it has to be completed ACROSS + RANKS: each boundary facet is labelled on exactly one rank, and a node on a + partition seam has its adjacent facets split between ranks (both edges are in the + local mesh, only one carries the label). Accumulating rank-locally would give a + seam node the wrong normal and make the answer rank-count dependent. The weighted + contributions are therefore summed over the DMPlex point SF 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. + + The analytic-normal path is NOT reduced: it evaluates a function of the node + coordinate, so every rank already computes the same value and summing copies would + only rescale it (and by a rank-count-dependent factor). """ dm = solver.dm dim = solver.mesh.dim @@ -119,15 +176,15 @@ def _boundary_velocity_nodes(solver, boundary, normal=None): cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) v0, v1 = dm.getDepthStratum(0) lsec = dm.getLocalSection() - interior_ref = cvec.mean(axis=0) # Boundary facets via the consolidated "UW_Boundaries" label (per-boundary labels do # not survive mesh adaptation); raises a clear error for an unknown boundary name. # In parallel a rank may own NO part of this boundary → a null IS; guard and return # no local nodes (calling getIndices() on a null IS would segfault). + # A rank owning no part of this boundary gets a null/empty IS and contributes + # nothing — but it must NOT return early: the cross-rank sum below is collective. sis = _boundary_stratum_is(dm, solver.mesh, boundary) - if not (sis and sis.getSize() > 0): - return [] - facets = [int(z) for z in sis.getIndices()] + facets = ([int(z) for z in sis.getIndices()] + if (sis and sis.getSize() > 0) else []) fS, fE = dm.getHeightStratum(1) # facets (edges in 2D, faces in 3D) def coord(q): @@ -160,18 +217,49 @@ def coord(q): else: const_normal = np.asarray(normal, dtype=float).ravel() - nacc = {} - pts = set() + nacc = {} # velocity node → Σ_f measure · n̂_f for f in facets: if not (fS <= f < fE): continue - # facet outward normal + # 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) + wgt = 1.0 if normal is None: - _, cent, nrm = dm.computeCellGeometryFVM(f) + vol, cent, nrm = dm.computeCellGeometryFVM(f) ne = np.asarray(nrm, dtype=float) ne = ne / (np.linalg.norm(ne) + 1e-30) - if np.dot(ne, np.asarray(cent) - interior_ref) < 0: - ne = -ne + # 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) # 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): @@ -184,12 +272,101 @@ def coord(q): else: ne = const_normal.copy() ne = ne / (np.linalg.norm(ne) + 1e-30) - nacc[q] = nacc.get(q, np.zeros(dim)) + ne - pts.add(q) - out = [] - for q in pts: - nrm = nacc[q] / (np.linalg.norm(nacc[q]) + 1e-30) - out.append((q, nrm)) + nacc[q] = nacc.get(q, np.zeros(dim)) + wgt * ne + + if normal is None: + nacc = _sum_facet_normals_across_ranks(solver, nacc) + + return [(q, v / (np.linalg.norm(v) + 1e-30)) + for q, v in sorted(nacc.items()) + if np.linalg.norm(v) > 0.0] + + +def _sum_facet_normals_across_ranks(solver, contribs): + """Complete the per-node facet sum across ranks: ``{point: Σ_f |f| n̂_f}`` restricted + to this rank's facets in, the sum over EVERY rank's facets out — the same value on + every rank that holds the node. + + Each boundary facet is labelled on exactly one rank, and on the rank that owns it + (measured: no facet is labelled twice and none is labelled away from its owner), so + a plain ADD is exact and needs no de-duplication. The velocity field already has + ``dim`` DOFs at exactly these points, so the sum rides the DM's own local↔global + scatter: ADD into the global vector accumulates the ghost copies onto the owner, + and scattering back gives every rank the identical total. + + COLLECTIVE — a rank owning no facet of this boundary still takes part. + """ + dm = solver.dm + dim = solver.mesh.dim + if dm.comm.getSize() == 1: + return contribs + + lsec = dm.getLocalSection() + lvec = dm.getLocalVec() + gvec = dm.getGlobalVec() + try: + lvec.set(0.0) + larr = lvec.getArray() + for q, v in contribs.items(): + lo = lsec.getFieldOffset(q, _VELOCITY_FIELD) + larr[lo:lo + dim] = v + gvec.set(0.0) + dm.localToGlobal(lvec, gvec, addv=PETSc.InsertMode.ADD_VALUES) + dm.globalToLocal(gvec, lvec) + summed = lvec.getArray() + out = {} + for q in _local_boundary_candidates(dm, lsec) | set(contribs): + lo = lsec.getFieldOffset(q, _VELOCITY_FIELD) + w = np.array(summed[lo:lo + dim], dtype=float) + # LOAD-BEARING: the candidate set over-collects (see + # _local_boundary_candidates); a point off this boundary contributes + # nothing anywhere, so its summed normal is exactly zero and it is + # dropped here. This is what makes the over-collection safe. + if w.any(): + out[q] = w + elif q in contribs: + # Velocity DOFs constrained out of the global vector (an essential + # BC meeting this boundary) come back zero. Keep this rank's own + # contribution so the node set is exactly what it always was. + # build_rotation drops these nodes (l2g < 0); boundary_normal_traction + # does NOT — it reads n̂·r_c at every returned node — so on such a node + # it reports the reaction of the ESSENTIAL constraint along a + # rank-locally accumulated normal. That was true before this change + # too, and the nodes are the corners where a rotated wall meets an + # essential one (measured: 2 of 75 on the test_1070 configuration, at + # every rank count, and never shared between ranks there). + out[q] = contribs[q] + return out + finally: + dm.restoreLocalVec(lvec) + dm.restoreGlobalVec(gvec) + + +def _local_boundary_candidates(dm, lsec): + """Velocity points on a support-1 facet of this rank's LOCAL mesh — a SUPERSET of + its boundary nodes, deliberately. + + The labelled facet list is not enough to enumerate a rank's boundary nodes: a rank + can OWN a node every one of whose labelled facets lives on a neighbour (the label + goes to one rank per facet, so a seam node's facets are split between them), and + such a node would otherwise never get a constraint row. + + Support 1 does NOT mean "on the domain boundary": the facets on the outer fringe of + the overlap layer also have support 1 locally while being interior globally, and + measured they are 13-27% of what this returns. The caller is what makes the result + correct — those points receive no contribution from any labelled facet, so their + summed normal is exactly zero and the ``w.any()`` test in + :func:`_sum_facet_normals_across_ranks` drops them. That test is load-bearing, not + an aside; do not remove it while over-collecting here. + """ + fS, fE = dm.getHeightStratum(1) + out = set() + for f in range(fS, fE): + if dm.getSupportSize(f) != 1: + continue + for q in (int(c) for c in dm.getTransitiveClosure(f)[0]): + if lsec.getFieldDof(q, _VELOCITY_FIELD) > 0: + out.add(q) return out @@ -1829,6 +2006,15 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="auto"): σ_nn is recovered from the CARTESIAN nodal reaction r_c = A·u − b: the nodal load is R_i = n̂_i · r_c(node_i), where n̂_i is THIS boundary's outward normal at node i. + + SIGN. ``n̂_i`` is whatever :func:`_boundary_velocity_nodes` produced for this + boundary, so σ_nn is measured along that direction. The geometric default is the + DOMAIN's outward normal everywhere, which on a concave boundary (an annulus or + shell inner arc) points toward the centre of curvature — the opposite of what UW3 + produced before #560, so σ_nn and the dynamic topography built on it reverse sign + there against earlier releases. An analytic ``normal=`` is used exactly as given + and is not reoriented, so ``X/|X|`` on an inner arc yields the other sign. See + "Which way is outward" in ``docs/developer/subsystems/rotated-freeslip.md``. Projecting the Cartesian reaction onto the boundary normal (rather than reading the rotated frame's normal row) is corner-correct — at a node shared with another rotated-free-slip boundary the rotated frame's first row is a mix of both walls' diff --git a/tests/parallel/test_1063_constrained_freeslip_parallel.py b/tests/parallel/test_1063_constrained_freeslip_parallel.py index 745f45b5..843f99a0 100644 --- a/tests/parallel/test_1063_constrained_freeslip_parallel.py +++ b/tests/parallel/test_1063_constrained_freeslip_parallel.py @@ -121,6 +121,15 @@ def _solve_gauge_diagnostics(): 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) 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 @@ -146,6 +155,17 @@ def test_constrained_raw_gauge_partition_independent(): 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) @pytest.mark.parametrize("kind", ["iso", "ti"]) def test_constrained_freeslip_partition_independent(kind): """The parallel solve must reproduce the serial reference: velocity bit- diff --git a/tests/parallel/test_1064_rotated_freeslip_parallel.py b/tests/parallel/test_1064_rotated_freeslip_parallel.py index 1ebe04ed..63a935a1 100644 --- a/tests/parallel/test_1064_rotated_freeslip_parallel.py +++ b/tests/parallel/test_1064_rotated_freeslip_parallel.py @@ -402,6 +402,18 @@ def test_rotated_freeslip_box_partition_independent(): 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) 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.""" @@ -454,6 +466,18 @@ def test_rotated_freeslip_mesh_owned_fmg_pickup(): 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) 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 @@ -468,6 +492,18 @@ def test_rotated_freeslip_spherical3d_partition_independent(): 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) def test_rotated_freeslip_spherical3d_topography_partition_independent(): """3D boundary-mass recovery gives partition-independent topography coefficients.""" coefficients = _spherical3d_topography_diagnostics() diff --git a/tests/parallel/test_1066_rotated_datum_parallel.py b/tests/parallel/test_1066_rotated_datum_parallel.py index f15e9ef8..65242134 100644 --- a/tests/parallel/test_1066_rotated_datum_parallel.py +++ b/tests/parallel/test_1066_rotated_datum_parallel.py @@ -27,6 +27,19 @@ _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(): RI, RO = 0.5, 1.0 mesh = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.1, qdegree=3) diff --git a/tests/parallel/test_1068_rotated_nodal_normal_parallel.py b/tests/parallel/test_1068_rotated_nodal_normal_parallel.py new file mode 100644 index 00000000..4244a036 --- /dev/null +++ b/tests/parallel/test_1068_rotated_nodal_normal_parallel.py @@ -0,0 +1,133 @@ +"""Parallel regression for the measure-weighted rotated nodal normal (issue #560). + +The nodal normal is a sum over ALL facets meeting the node, so it only survives +partitioning if the sum is completed across ranks. Each boundary facet is labelled on +exactly one rank, so a node on a partition seam sees only a SUBSET of its facets +locally: accumulating rank-locally gives a seam node a different normal on each rank +and makes the answer rank-count dependent. Two further traps live here — a rank can OWN +a boundary node whose labelled facets are all on neighbours, and an outward-orientation +test taken against the mean of the rank's own coordinates can orient two facets of the +same node oppositely so that they CANCEL in the sum. + +Measured on the skewed annulus (both arcs exactly circular, nodes moved to non-uniform +angles so only the facet LENGTHS change), geometric normal, ``|A z|/|A|_F``: + + np before after + 1 1.541e-07 6.624e-20 + 2 1.611e-06 6.623e-20 + 4 1.392e-05 6.690e-20 + +Run with:: + + mpirun -n 2 python -m pytest --with-mpi tests/parallel/test_1068_rotated_nodal_normal_parallel.py + mpirun -n 4 python -m pytest --with-mpi tests/parallel/test_1068_rotated_nodal_normal_parallel.py +""" +import numpy as np +import pytest +import sympy +import underworld3 as uw +import underworld3.utilities.rotated_bc as rbc + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.timeout(300)] + + +def _skewed_annulus(tag, skew=0.04): + mesh = uw.meshing.Annulus(radiusOuter=1.0, radiusInner=0.55, cellSize=0.15, + qdegree=3) + c = mesh.X.coords.copy() + th = np.arctan2(c[:, 1], c[:, 0]) + r = np.hypot(c[:, 0], c[:, 1]) + th = th + skew * np.sin(3.0 * th) + mesh.deform(np.column_stack([r * np.cos(th), r * np.sin(th)])) + v = uw.discretisation.MeshVariable(f"vN{tag}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"pN{tag}", mesh, 1, degree=1, continuous=False) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + x, y = mesh.X + s.bodyforce = sympy.Matrix([[x, y]]) * sympy.cos(4 * sympy.atan2(y, x)) + s.penalty = 0.0 + s.tolerance = 1e-9 + s.add_rotated_freeslip_bc(0, "Upper") # geometric normal + s.add_rotated_freeslip_bc(0, "Lower") + s.petsc_use_pressure_nullspace = True + return s + + +class _ConstantPressureResidual: + """``|A z|/|A|_F`` for the attached constant pressure, by PETSc matvec only. + + The Frobenius norm is not sigma_max, so this is NOT comparable with the serial + dense-SVD number in ``tests/test_1018_rotated_nodal_normal.py`` — it is only + comparable with itself across rank counts, which is exactly what is at issue.""" + + def __enter__(self): + self.value = None + self._orig = rbc._solve_rotated_iterative + + def _capture(solver, Ahat, bhat, Q, Qt, normal_rows, **kw): + if self.value is None and kw.get("nsp") is not None: + z = kw["nsp"].getVecs()[0].copy() + z.normalize() + r = Ahat.createVecLeft() + Ahat.mult(z, r) + self.value = r.norm() / Ahat.norm(2) # 2 == NORM_FROBENIUS + z.destroy() + r.destroy() + return self._orig(solver, Ahat, bhat, Q, Qt, normal_rows, **kw) + + rbc._solve_rotated_iterative = _capture + return self + + def __exit__(self, *exc): + rbc._solve_rotated_iterative = self._orig + return False + + +def test_nodal_normal_is_bit_identical_across_ranks(): + """Every rank holding a shared boundary node must compute the SAME normal, to the + last bit — otherwise the constraint frame depends on the partition.""" + s = _skewed_annulus("B") + s.solve() + dm, dim = s.dm, s.mesh.dim + csec = dm.getCoordinateSection() + cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) + v0, v1 = dm.getDepthStratum(0) + + local = [] + for wall in ("Upper", "Lower"): + for q, nrm in rbc._boundary_velocity_nodes(s, wall): + c = rbc._point_coord(dm, dim, cvec, csec, v0, v1, q) + local.append((wall, tuple(round(float(t), 9) for t in c), + tuple(float(t) for t in nrm))) + + gathered = uw.mpi.comm.gather(local, root=0) + if uw.mpi.rank == 0: + seen = {} + shared = disagreeing = 0 + worst = 0.0 + for rows in gathered: + for wall, coord, nrm in rows: + key = (wall, coord) + if key in seen: + shared += 1 + d = max(abs(a - b) for a, b in zip(seen[key], nrm)) + worst = max(worst, d) + disagreeing += d > 0.0 + else: + seen[key] = nrm + assert shared > 0, "no boundary node is shared — the test proves nothing" + assert disagreeing == 0, ( + f"{disagreeing} of {shared} shared boundary nodes disagree between " + f"ranks, max |dn| = {worst:.3e}; the nodal normal is partition dependent") + + +def test_constant_pressure_is_a_null_vector_in_parallel(): + """The serial result must survive partitioning: rank-locally accumulated normals + read 1.6e-06 at np=2 and 1.4e-05 at np=4 against 6.6e-20 here.""" + with _ConstantPressureResidual() as cap: + _skewed_annulus("R").solve() + assert cap.value is not None, "the operator capture never fired" + assert cap.value < 1e-15, ( + f"|A z|/|A|_F = {cap.value:.3e} at np={uw.mpi.size}; the constant pressure is " + "not a null vector of the constrained operator (#560)") diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index 6aa789ab..53af6ea4 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -249,23 +249,6 @@ def __str__(self): f"radius_changed={self.radius_changed}") -def _deformed_box_rotated_solve(tag, coord_scale=1.0): - """A rotated free-slip solve on the deformed box, optionally with every - coordinate multiplied by ``coord_scale``. Used to MEASURE the direction the - system does not determine (#560): the physical response to a coordinate - change of a few machine epsilons is of that order, so anything the solution - does beyond that is the unpinned mode.""" - mesh = uw.meshing.StructuredQuadBox( - elementRes=(10, 10), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3) - c = mesh.X.coords.copy() - c[:, 1] += 0.02 * c[:, 1] * np.sin(np.pi * c[:, 0]) - mesh.deform(c * coord_scale) - k = uw.function.expression(r"k_probe", 1.0, "viscosity") - s, v = _rampable_rotated_stokes(mesh, k, tag) - s.solve() - return np.asarray(v.data).copy() - - def _rigid_body_decomposition(coords, diff): """Split a velocity difference field into rigid-body modes and the rest. @@ -411,24 +394,19 @@ def test_rotated_workspace_deform_invalidates(): the deform, and both solves converged. This is the contract #543 wrote the test for and it is asserted hard. - **(b) the two answers agree.** Rotated free-slip on a CURVED boundary loses - the constant-pressure gauge (#560), leaving one direction the operator does - not pin: a coordinate change of two machine epsilons moves the velocity by - 1.3e-01, and the move does not scale with the perturbation. So "the two - solves agree" is false as stated — a plain comparison passes only where the - two assemblies happen to agree bitwise, which is why this test was green on - macOS (err exactly 0.0 in 81 consecutive runs across two PETSc toolchains - and nine PYTHONHASHSEEDs) and intermittently red on CI (err 7e-2 to 1.2e-1, - the size of the unpinned component rather than a drift). - - The unpinned subspace is one-dimensional (five perturbations move the - answer along the same direction to cosine 1.000000), so the claim is - narrowed rather than dropped: measure that direction with one extra - perturbed solve and require agreement in every other direction. The - tolerance is NOT relaxed — 1e-6, as before — it is the claim that is made - honest. When #560 is fixed the probe stops finding a direction, the - projection becomes a no-op, and the test compares the solutions directly - again with no further edit. + **(b) the two answers agree.** This claim WAS ill posed. Rotated free-slip + on a CURVED boundary used to lose the constant-pressure gauge (#560), + leaving one direction the operator did not pin: a coordinate change of two + machine epsilons moved the velocity by 1.3e-01 and the move did not scale + with the perturbation, so a plain comparison passed only where the two + assemblies happened to agree bitwise — green on macOS (err exactly 0.0 in + 81 consecutive runs) and intermittently red on CI (err 7e-2 to 1.2e-1, the + size of the unpinned component rather than a drift). The test then compared + the two solutions with that one direction projected out. + + #560 is fixed and the same probe now measures 9.9e-15, so the projection + and the probe solve are gone and the two solutions are compared directly, + at the same 1e-6 tolerance as before. The instrumentation below (gauge decisions, constrained-row counts, locator tallies, rigid-body decomposition) stays: it is what turned an unreadable @@ -531,75 +509,35 @@ def test_rotated_workspace_deform_invalidates(): # ------------------------------------------------------------------ # (b) the answers agree, in the directions the system actually determines # ------------------------------------------------------------------ - # Rotated free-slip on a CURVED boundary loses the constant-pressure gauge - # (#560): the pressure level runs to ~1e4 against a pressure variation of - # 6.6e-2, and the solution acquires a component along one unpinned - # direction whose amplitude is set by round-off. Measured: a coordinate - # change of two machine epsilons (4.44e-16) moves the velocity by 1.33e-01, - # and the size of the move does not track the size of the perturbation - # (4.4e-16, 2.2e-15, 1e-14, 1e-12 and 1e-9 all give 4e-2 to 2e-1). The same - # solve on a STRAIGHT-walled box moves by 5.9e-15, and native essential - # free-slip on this same deformed mesh moves by 4.8e-13 — so it is the - # rotated path on a curved boundary, and it is present at the merge base. + # This comparison used to be ill posed, and is not any more. Rotated + # free-slip on a CURVED boundary lost the constant-pressure gauge (#560), + # leaving one direction the operator did not pin: a coordinate change of + # two machine epsilons moved the velocity by 1.33e-01 and the size of the + # move did not track the size of the perturbation, so "the two solves + # agree" was false as stated. The test passed only where the two + # assemblies happened to agree bitwise — green on macOS/arm64, red on CI + # at err 7e-2 to 1.2e-1, the SIZE of the unpinned component rather than a + # drift. It compared the two solutions with the unpinned direction + # projected out, measured by an extra perturbed solve. # - # So "the two solves agree" is not a property this system has, and a plain - # comparison passes only where the two assemblies happen to agree bitwise - # (macOS/arm64: err is exactly 0.0 in 81 consecutive runs across two PETSc - # toolchains and nine PYTHONHASHSEEDs; CI's Linux build: err ~7e-2 to - # 1.2e-1, which is the SIZE of the unpinned component, not a drift). - # - # The unpinned subspace is exactly ONE-dimensional — five different - # perturbations move the answer along the same direction to cosine - # 1.000000, the normalised difference set has singular values - # [2.236, 4.4e-9, ...], and removing the leading direction leaves 2e-9 of - # each difference. So the comparison can be made well posed rather than - # abandoned: measure that direction with one extra perturbed solve and - # assert the two solutions agree in every OTHER direction. - probe = _deformed_box_rotated_solve("Prb", coord_scale=1.0 + 2 * _EPS) - unpinned = np.asarray(probe) - np.asarray(v_c.data) - unpinned_size = np.linalg.norm(unpinned) / np.linalg.norm(v_c.data) - - if unpinned_size > 1.0e-3: - # #560 is present (the expected branch today). Project it out. - direction = (unpinned / np.linalg.norm(unpinned)).ravel() - flat = diff.ravel() - residual = flat - float(np.dot(flat, direction)) * direction - constrained_err = np.linalg.norm(residual) / np.linalg.norm(v_c.data) - branch = (f"#560 present: the unpinned direction carries " - f"{unpinned_size:.3e} of the solution, projected out") - else: - # #560 has been fixed — there is no unpinned direction to remove, so - # compare directly and let this test go back to its full strength. - constrained_err = err - branch = (f"#560 appears FIXED (a 2-eps perturbation moves the answer " - f"by only {unpinned_size:.3e}) — the projection below is " - f"now a no-op and this test is comparing solutions directly. " - f"Delete the projection and the _deformed_box_rotated_solve " - f"probe.") - - # NEGATIVE CONTROL: the projection must not absorb a genuine discrepancy. - # Inject a difference orthogonal to the unpinned direction and check it - # survives, or "constrained_err is small" would be true of anything. + # #560 is fixed (the nodal normal is measure-weighted, so the constant + # pressure is a null vector of the constrained operator again) and that + # probe now measures 9.9e-15 rather than 1.3e-01. The projection and the + # probe solve are therefore gone and the solutions are compared directly, + # at the original 1e-6 tolerance. + constrained_err = err + + # NEGATIVE CONTROL: inject a discrepancy and check the comparison sees it, + # or "constrained_err is small" would be true of anything. injected = np.random.default_rng(6).normal(size=diff.shape).ravel() - if unpinned_size > 1.0e-3: - d_hat = (unpinned / np.linalg.norm(unpinned)).ravel() - injected -= float(np.dot(injected, d_hat)) * d_hat injected *= 1.0e-3 * np.linalg.norm(v_c.data) / np.linalg.norm(injected) - poisoned = diff.ravel() + injected - if unpinned_size > 1.0e-3: - d_hat = (unpinned / np.linalg.norm(unpinned)).ravel() - poisoned = poisoned - float(np.dot(poisoned, d_hat)) * d_hat - poisoned_err = np.linalg.norm(poisoned) / np.linalg.norm(v_c.data) + poisoned_err = np.linalg.norm(diff.ravel() + injected) / np.linalg.norm(v_c.data) assert poisoned_err > 1.0e-4, ( - f"a deliberate 1e-3 discrepancy orthogonal to the unpinned direction " - f"survives the projection as only {poisoned_err:.3e}, so the " - f"constrained comparison below would not notice a real disagreement") + f"a deliberate 1e-3 discrepancy shows up as only {poisoned_err:.3e}, so " + f"the comparison below would not notice a real disagreement") assert constrained_err < 1e-6, ( - f"post-deform solve differs from fresh control by {constrained_err:.2e} " - f"OUTSIDE the direction the system leaves undetermined " - f"(raw difference {err:.2e})\n" - f" branch : {branch}\n" + f"post-deform solve differs from fresh control by {constrained_err:.2e}\n" f" first solve : {report_1}\n" f" post-deform : {report_deformed}\n" f" fresh control : {report_control}\n" diff --git a/tests/test_1018_rotated_nodal_normal.py b/tests/test_1018_rotated_nodal_normal.py new file mode 100644 index 00000000..1eda1570 --- /dev/null +++ b/tests/test_1018_rotated_nodal_normal.py @@ -0,0 +1,343 @@ +"""The rotated free-slip nodal normal must be consistent with the assembly (issue #560). + +A boundary node on more than one facet gets ONE nodal normal while the assembler +integrates the boundary term facet by facet. If that normal is the plain bisector of +its facet normals, a kinked boundary with unequal facets leaves a residual in the +node's FREE tangential row, the exact constant-pressure vector stops being a null +vector of the constrained rotated operator, and the pressure gauge goes unpinned. + +The measurable is ``|A z| / |A|`` with ``z`` the constant-pressure vector UW3 attaches +as the null space and ``A`` the constrained operator handed to the KSP. It is the right +metric because it scales smoothly with the deformation (amplitude cubed) — ``mean(p)`` +and the run-to-run velocity move are the ROUND-OFF amplitude of an unpinned direction +and are chaotic, so they report the presence of the defect but not its size. + +Each check carries its negative control: the same measurement with the pre-fix +bisector normal, which must FAIL the tolerance. Without that the tests would pass on +a metric that cannot see the defect. +""" +import numpy as np +import pytest +import sympy +import underworld3 as uw +import underworld3.utilities.rotated_bc as rbc + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + +_MACHINE = 1e-15 # |A z|/|A| at machine level; the straight box sits at 3e-18 + + +def _bisector_nodes(solver, boundary, normal=None, + _real=rbc._boundary_velocity_nodes): + """The pre-#560 accumulation, frozen here as the negative control. + + Identical to ``rbc._boundary_velocity_nodes`` except that each facet normal is + normalised to unit length BEFORE accumulating, so a node's normal is the bisector + of its facet normals rather than the measure-weighted average. The analytic-normal + path is delegated to the real function — it is unchanged by the fix. + + ``_real`` binds the genuine function AT IMPORT. Looking it up at call time would + self-recurse whenever this control is installed over the module attribute and an + analytic normal is passed. + + It also freezes the OLD outward test (against the mean of the mesh coordinates), + so on a CONCAVE geometric boundary it differs from the real function by a SIGN as + well as by the weighting — which is what + ``test_geometric_normal_points_out_of_the_domain_on_a_concave_boundary`` uses it + for. Do not extend a bit-identity comparison onto an inner arc's geometric normal + expecting agreement: the sign difference there is the intended change, not a + regression. + """ + if normal is not None: + return _real(solver, boundary, normal=normal) + dm = solver.dm + dim = solver.mesh.dim + cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) + lsec = dm.getLocalSection() + interior_ref = cvec.mean(axis=0) + sis = rbc._boundary_stratum_is(dm, solver.mesh, boundary) + if not (sis and sis.getSize() > 0): + return [] + fS, fE = dm.getHeightStratum(1) + nacc, pts = {}, set() + for f in (int(z) for z in sis.getIndices()): + if not (fS <= f < fE): + continue + _, cent, nrm = dm.computeCellGeometryFVM(f) + ne = np.asarray(nrm, dtype=float) + ne = ne / (np.linalg.norm(ne) + 1e-30) + if np.dot(ne, np.asarray(cent) - interior_ref) < 0: + ne = -ne + for q in (int(c) for c in dm.getTransitiveClosure(f)[0]): + if lsec.getFieldDof(q, rbc._VELOCITY_FIELD) <= 0: + continue + nacc[q] = nacc.get(q, np.zeros(dim)) + ne + pts.add(q) + return [(q, nacc[q] / (np.linalg.norm(nacc[q]) + 1e-30)) for q in pts] + + +class _CaptureOperator: + """Capture the constrained rotated operator and its attached null space. + + The operator only exists inside the rotated solve, so the linear-solve entry point + is wrapped for the duration of the ``with`` block. Nothing else can see it. + """ + + def __init__(self): + self.A = None + self.nsp = None + + def __enter__(self): + self._orig = rbc._solve_rotated_iterative + + def _capture(solver, Ahat, bhat, Q, Qt, normal_rows, **kw): + if self.A is None: + n = Ahat.getSize()[0] + ai, aj, av = Ahat.getValuesCSR() # never MatConvert: it breaks the PC + dense = np.zeros((n, n)) + for r in range(n): + dense[r, aj[ai[r]:ai[r + 1]]] = av[ai[r]:ai[r + 1]] + self.A = dense + self.nsp = [np.asarray(v.array_r).copy() + for v in kw["nsp"].getVecs()] + return self._orig(solver, Ahat, bhat, Q, Qt, normal_rows, **kw) + + rbc._solve_rotated_iterative = _capture + return self + + def __exit__(self, *exc): + rbc._solve_rotated_iterative = self._orig + return False + + @property + def constant_pressure_residual(self): + """``|A z| / |A|`` for the attached constant-pressure vector.""" + z = self.nsp[0] / np.linalg.norm(self.nsp[0]) + return np.linalg.norm(self.A @ z) / np.linalg.norm(self.A, 2) + + +class _BisectorNormals: + """Run the enclosed solves with the pre-#560 bisector nodal normal.""" + + def __enter__(self): + self._orig = rbc._boundary_velocity_nodes + rbc._boundary_velocity_nodes = _bisector_nodes + return self + + def __exit__(self, *exc): + rbc._boundary_velocity_nodes = self._orig + return False + + +def _deformed_box(tag, amplitude, res=10): + """Configuration B of issue #560: a box whose top is bowed by + ``amplitude * y * sin(pi x)`` — a kinked rotated wall with unequal facets.""" + mesh = uw.meshing.StructuredQuadBox( + elementRes=(res, res), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3) + if amplitude != 0.0: + c = mesh.X.coords.copy() + c[:, 1] += amplitude * c[:, 1] * np.sin(np.pi * c[:, 0]) + mesh.deform(c) + return _rotated_box_solver(tag, mesh) + + +def _rotated_box_solver(tag, mesh): + v = uw.discretisation.MeshVariable(f"vK{tag}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"pK{tag}", mesh, 1, degree=1, continuous=False) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + x, y = mesh.X + s.bodyforce = sympy.Matrix([[0.0, sympy.sin(sympy.pi * x) * sympy.cos(sympy.pi * y)]]) + s.penalty = 0.0 + s.tolerance = 1e-9 + for wall in ("Top", "Bottom", "Left", "Right"): + s.add_rotated_freeslip_bc(0, wall) + s.petsc_use_pressure_nullspace = True + return s + + +def _skewed_annulus(tag, skew): + """Both boundaries stay exactly circular; the nodes are moved to non-uniform + angles, so only the facet LENGTHS change. Geometric (default) normals.""" + mesh = uw.meshing.Annulus(radiusOuter=1.0, radiusInner=0.55, cellSize=0.15, + qdegree=3) + if skew != 0.0: + c = mesh.X.coords.copy() + th = np.arctan2(c[:, 1], c[:, 0]) + r = np.hypot(c[:, 0], c[:, 1]) + th = th + skew * np.sin(3.0 * th) + mesh.deform(np.column_stack([r * np.cos(th), r * np.sin(th)])) + v = uw.discretisation.MeshVariable(f"vA{tag}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"pA{tag}", mesh, 1, degree=1, continuous=False) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + x, y = mesh.X + s.bodyforce = sympy.Matrix([[x, y]]) * sympy.cos(4 * sympy.atan2(y, x)) + s.penalty = 0.0 + s.tolerance = 1e-9 + s.add_rotated_freeslip_bc(0, "Upper") + s.add_rotated_freeslip_bc(0, "Lower") + s.petsc_use_pressure_nullspace = True + return s + + +def _residual(build, *args): + with _CaptureOperator() as cap: + build(*args).solve() + return cap.constant_pressure_residual + + +def test_deformed_box_constant_pressure_stays_a_null_vector(): + """Configuration B: the attached constant pressure must annihilate the constrained + operator on a kinked rotated wall, as it already does on a straight one.""" + straight = _residual(_deformed_box, "S", 0.0) + deformed = _residual(_deformed_box, "D", 0.02) + assert straight < _MACHINE, f"straight box control |Az|/|A| = {straight:.3e}" + assert deformed < _MACHINE, ( + f"deformed rotated free-slip box leaves |Az|/|A| = {deformed:.3e}; the " + "constant pressure is not a null vector and the gauge is unpinned (#560)") + + +def test_bisector_normal_fails_the_same_measurement(): + """Negative control: with the pre-fix bisector normal the metric must FIRE. + + Without this the test above could pass on a measurement that cannot see the + defect at all.""" + with _BisectorNormals(): + deformed = _residual(_deformed_box, "N", 0.02) + straight = _residual(_deformed_box, "M", 0.0) + assert deformed > 1e-11, ( + f"the bisector normal gave |Az|/|A| = {deformed:.3e}; the metric is not " + "detecting the defect it is supposed to detect") + assert straight < _MACHINE, ( + "the straight box is clean under BOTH normals — the negative control must " + "isolate the kink, not the whole rotated path") + + +def test_skewed_annulus_constant_pressure_stays_a_null_vector(): + """A curved rotated wall is clean while its facets are equal and acquires the + defect as soon as they are not — the case that separates the mechanism (unequal + kinked facets) from the symptom (a deformed box).""" + skewed = _residual(_skewed_annulus, "K", 0.04) + with _BisectorNormals(): + skewed_bisector = _residual(_skewed_annulus, "B", 0.04) + assert skewed < _MACHINE, ( + f"skewed annulus leaves |Az|/|A| = {skewed:.3e} (#560)") + assert skewed_bisector > 1e-11, ( + f"negative control did not fire on the skewed annulus " + f"({skewed_bisector:.3e})") + + +def test_flat_walls_and_analytic_normals_are_bit_identical(): + """An AXIS-ALIGNED wall and the analytic-normal override must be unchanged to the + last bit: the first has facet normals with exactly 0/±1 components, so the weights + cancel in the normalisation whatever they are; the second never consults a facet. + (A flat but TILTED wall is not covered — there the weighting moves ~1 ulp.)""" + straight = _deformed_box("I", 0.0) + straight.solve() + for wall in ("Top", "Bottom", "Left", "Right"): + fixed = dict(rbc._boundary_velocity_nodes(straight, wall)) + old = dict(_bisector_nodes(straight, wall)) + assert fixed and fixed.keys() == old.keys() + for q in fixed: + assert np.array_equal(fixed[q], old[q]), ( + f"straight wall {wall!r} node {q} normal moved: " + f"{fixed[q]} vs {old[q]}") + + annulus = _skewed_annulus("J", 0.04) + annulus.solve() + x, y = annulus.mesh.X + radial = sympy.Matrix([[x, y]]) / sympy.sqrt(x**2 + y**2) + for arc in ("Upper", "Lower"): + fixed = dict(rbc._boundary_velocity_nodes(annulus, arc, normal=radial)) + old = dict(_bisector_nodes(annulus, arc, normal=radial)) + # without this guard an empty stratum would pass with zero comparisons, + # and this is the only test exercising the analytic path at all + assert fixed and fixed.keys() == old.keys() + for q in fixed: + assert np.array_equal(fixed[q], old[q]), ( + f"analytic normal on {arc!r} node {q} moved") + + # and the control: on the KINKED geometric path the two must differ, or the + # comparison above is vacuous. + kinked = _deformed_box("H", 0.02) + kinked.solve() + fixed = dict(rbc._boundary_velocity_nodes(kinked, "Top")) + old = dict(_bisector_nodes(kinked, "Top")) + assert any(not np.array_equal(fixed[q], old[q]) for q in fixed), ( + "the deformed top gave identical normals under both rules — the " + "bit-identity checks above prove nothing") + + +def test_geometric_normal_points_out_of_the_domain_on_a_concave_boundary(): + """The geometric normal is the DOMAIN's outward normal, which on a concave + boundary points toward the centre of curvature. + + This pins a user-visible convention. Before #560 the outward test was "away from + the mean of the mesh coordinates", which on an annulus inner arc points INTO the + domain; it is now "away from the facet's own support cell", which is outward + everywhere. σ_nn and dynamic topography read on an inner boundary through the + geometric normal therefore changed sign, and without this assertion nothing would + stop that flipping back silently. + + An analytic ``normal=`` is used exactly as supplied and is NOT reoriented, so + ``X/|X|`` on an inner arc is inward-of-domain — asserted here too, because the + disagreement between the two paths is deliberate and is what the docs describe.""" + annulus = _skewed_annulus("C", 0.0) + annulus.solve() + dm, dim = annulus.dm, annulus.mesh.dim + csec = dm.getCoordinateSection() + cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) + v0, v1 = dm.getDepthStratum(0) + + def radial_components(boundary, normal=None, nodes=None): + out = [] + source = nodes or rbc._boundary_velocity_nodes + for q, nrm in source(annulus, boundary, normal=normal): + x = np.asarray(rbc._point_coord(dm, dim, cvec, csec, v0, v1, q)) + out.append(float(np.dot(nrm, x / np.linalg.norm(x)))) + return np.array(out) + + outer = radial_components("Upper") + inner = radial_components("Lower") + assert outer.size and inner.size + assert np.all(outer > 0.99), ( + f"outer arc normal is not outward (min radial component {outer.min():.6f})") + assert np.all(inner < -0.99), ( + f"inner arc normal must point OUT OF THE DOMAIN, i.e. toward the origin " + f"(max radial component {inner.max():.6f}); the pre-#560 rule gave +1 here") + + # NEGATIVE CONTROL: the frozen pre-#560 accumulation orients against the mean of + # the mesh coordinates, which on an inner arc is the WRONG way. If this did not + # come out +1 the assertion above would be pinning nothing. + old_inner = radial_components("Lower", nodes=_bisector_nodes) + assert np.all(old_inner > 0.99), ( + f"the pre-#560 rule should give an INWARD-of-domain normal on the inner arc " + f"(got min {old_inner.min():.6f}); without that the sign convention asserted " + f"above is not actually a change anything could regress to") + + x, y = annulus.mesh.X + radial = sympy.Matrix([[x, y]]) / sympy.sqrt(x**2 + y**2) + analytic_inner = radial_components("Lower", normal=radial) + assert np.all(analytic_inner > 0.99), ( + "an analytic X/|X| on the inner arc should be left INWARD-of-domain — the " + "override is applied exactly as given, and the docs say so") + + +@pytest.mark.level_2 +def test_residual_does_not_grow_with_deformation_amplitude(): + """Before the fix ``|A z|`` grew as amplitude cubed. After it, the deformation + must not drive it at all.""" + amplitudes = (0.0, 0.005, 0.02, 0.05) + fixed = [_residual(_deformed_box, f"A{k}", a) for k, a in enumerate(amplitudes)] + assert max(fixed) < _MACHINE, ( + "|Az|/|A| over the amplitude sweep " + + ", ".join(f"{a:g}: {r:.3e}" for a, r in zip(amplitudes, fixed))) + + with _BisectorNormals(): + control = [_residual(_deformed_box, f"C{k}", a) for k, a in enumerate(amplitudes)] + assert control[-1] > 100 * control[1], ( + "the negative control did not reproduce the growth with amplitude: " + + ", ".join(f"{a:g}: {r:.3e}" for a, r in zip(amplitudes, control)))