diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index 2a8f28f2..574bf172 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -6,6 +6,39 @@ This log tracks significant development work at a conceptual level, suitable for ## 2026 Q3 (July – September) +### The Free Surface Reaches the Spherical Shell (July 2026) + +**`uw.systems.FreeSurface` now runs in 3D on a spherical shell** — the same +exponential three-number integrator, held-lid σ_nn recovery and strong +material-boundary datum, with the surface machinery made dimension-general +rather than ported piecewise: + +- The datum gauge (mean removal) is an FE trace-mass reduction over + owned boundary facets — no ordered ring, no gather; the same code is the + 2D line gauge and the 3D area gauge. On the way it resolved a real 2D + defect: the deforming-ring strong-datum solves used to stall at a ~2e-3 + residual floor, which turned out to be three stacked causes (arc-length vs + FE trace weights; the datum's *directed* mean flux through the deformed + facet normals, now stripped with the same FE surface integral the residual + uses; and the constant-pressure gauge mode, which the inner solver projects + and the outer loop therefore now measures in the quotient space). With all + three closed, every step of the power-law acceptance run converges. +- σ_nn on a 3D P2 boundary is recovered by **P1 projection** (edge-midpoint + loads folded exactly onto vertices, sound P1 lumped triangle mass) — chosen + over the consistent P2 mass because its vertex-integral checkerboard sits + exactly at the vertices the P1 topography field consumes. +- The two genuinely 2D features (ring Taubin filter, tangential transport) + are refused explicitly in 3D; everything else is shared code. + +First 3D evidence (spherical Y20 topographic relaxation, constant-density +shell): exponential decay at an O(1) shell correction below the half-space +Cathles rate, in the physically correct direction, with the equilibrium +modal bias falling 16% → 2% of the initial amplitude over one resolution +step (the known discrete recovery defect, resolution-convergent). The +detailed benchmarking — analytic shell-rate comparison, convergence study, +low-Ra spherical convection, 3D parallel — is deliberately left to the +review pass. + ### One Rotated Free-Slip Path, Now With a Prescribed Wall-Normal Velocity (July 2026) **Rotated strong free-slip now takes a prescribed wall-normal velocity datum diff --git a/src/underworld3/systems/free_surface.py b/src/underworld3/systems/free_surface.py index 1fbbbb4d..700a1888 100644 --- a/src/underworld3/systems/free_surface.py +++ b/src/underworld3/systems/free_surface.py @@ -139,6 +139,13 @@ def __init__( self.composition = composition self._conserve_integrand = conserve self._smooth_length = smooth_length + # sigma_nn de-smear: "lumped" is the monotone 2D default; on a 3D P2 + # trace the lumped vertex mass is identically zero, and the consistent + # P2 mass carries the vertex-integral checkerboard (#404 hold) exactly + # at the vertices the P1 h_inf field reads — so 3D translates the + # default to the sound P1-PROJECTED recovery (boundary_flux mass="p1"). + if stokes.mesh.dim == 3 and mass == "lumped": + mass = "p1" self._mass = mass self.max_surface_cfl = max_surface_cfl # First-pass along-surface (tangential) transport of the surface fields: @@ -176,18 +183,29 @@ def __init__( "free surface the sigma_nn recovery zigzags. Prefer continuous pressure." ) - # Topography direction: vertical (last axis) on a Cartesian box, radial on a - # cylindrical annulus. The surface height and the mesh deformation follow it. - if self.mesh.dim != 2: - raise NotImplementedError( - "FreeSurface is 2D-only for now: the surface ring machinery (filter, " - "tangential transport, arc-length datum) and the sigma_nn de-smear have " - "no 3D counterparts yet (the 3D boundary mass is PR #404's scope). " - "Without this guard a 3D run would proceed and be silently wrong." - ) - self._cylindrical = ( - self.mesh.CoordinateSystem.coordinate_type == CoordinateSystemType.CYLINDRICAL2D - ) + # Topography direction: vertical (last axis) on a Cartesian box/slab, radial + # on a cylindrical annulus or spherical shell. The surface height and the + # mesh deformation follow it. `_radial` steers the geometry (any dimension); + # `_cylindrical` additionally selects the 2D ring's angular ordering. + ctype = self.mesh.CoordinateSystem.coordinate_type + self._cylindrical = ctype == CoordinateSystemType.CYLINDRICAL2D + self._radial = self._cylindrical or ctype == CoordinateSystemType.SPHERICAL + # 3D runs on the dimension-general primitives (facet trace-mass gauge, + # directed flux strip, unordered surface gather, nodal carrier). The two + # 2D-ONLY features are refused per piece rather than by a blanket guard: + if self.mesh.dim == 3: + if tangent_advect is not None: + raise NotImplementedError( + "FreeSurface: tangential surface transport is 2D-only (it runs " + "on the ordered surface ring; a 3D counterpart needs surface FE " + "advection). Construct with tangent_advect=None in 3D." + ) + if int(surface_filter) > 0: + raise NotImplementedError( + "FreeSurface: the Taubin surface filter is 2D-only (1-D ring " + "stencil); pass surface_filter=0 in 3D (the trace-mass gauge " + "and the interior carrier do not require it)." + ) self._walls = self._classify_walls() # Reference-configuration surface nodes: used ONCE to match each surface node to # its row in the mesh coordinate field and in the derived surface fields. Row @@ -252,15 +270,16 @@ def _ring_coords(self): def _surface_height(self, coords): """The coordinate along the topography direction: the last axis on a Cartesian - box, the radius on a cylindrical annulus.""" - if self._cylindrical: + box/slab, the radius on a cylindrical annulus or spherical shell.""" + if self._radial: return np.linalg.norm(coords, axis=1) return np.asarray(coords[:, -1], dtype=float) 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).""" - if self._cylindrical: + increment is deformed along — vertical (Cartesian) or radial (annulus / + spherical shell); dimension-general.""" + if self._radial: r = np.linalg.norm(coords, axis=1) r[r == 0.0] = 1.0 return coords / r[:, None] @@ -333,25 +352,73 @@ def _field_rows(self, field, coords): return rows[order], coords[order, 0] def _ring_weights(self): - r"""Arc-length quadrature weights on the globally s-sorted surface ring. - - The trapezoidal weight of node :math:`i` is half the distance to each neighbour - — :math:`\tfrac12 (s_{i+1}-s_{i-1}) r_i` on a cylindrical ring, :math:`\tfrac12 - (x_{i+1}-x_{i-1})` on an open Cartesian surface (with half-cells at the ends). - Radii are read live, so the weights follow the deforming surface. + r"""P1 lumped trace-mass quadrature weights on the globally s-sorted surface + ring: :math:`w_i = \oint \phi_i \,\mathrm{d}s`, accumulated per boundary FACET + from the DMPlex (each facet contributes measure/nverts to each of its + vertices — edge length/2 in 2D, triangle area/3 in 3D), from LIVE + coordinates so the weights follow the deforming surface. + + The facet (chord) measure — not an arc approximation — is deliberate: the + prescribed :math:`\tilde u_n` must be flux-free in the FINITE-ELEMENT sense + (:math:`\oint` over the deformed polygon the discretisation actually + integrates), or the strong datum asks the incompressible interior for a flow + that does not exist. With trapezoid arc weights the consistent solve floored + at rel ~2e-3 with the entire residual in the pressure (divergence) rows — + measured, block-split; the trace-mass weights are exact for the P1 datum + field by construction. Dimension-general: the same accumulation is the + area-weighted gauge on a 3D boundary triangulation. """ - s = self._s_sorted # UNIQUE ring: one entry per node - if s.size == 0: - return s - if self._ring_period is not None: - ds = 0.5 * np.mod(np.roll(s, -1) - np.roll(s, 1), self._ring_period) - radius = self._ring_gather(np.linalg.norm(self._ring_coords, axis=1)) - return ds * radius - ds = np.empty_like(s) - ds[1:-1] = 0.5 * (s[2:] - s[:-2]) - ds[0] = 0.5 * (s[1] - s[0]) - ds[-1] = 0.5 * (s[-1] - s[-2]) - return ds + return self._ring_gather(self._surface_weights_local(), op="sum") + + def _surface_weights_local(self): + """This rank's OWNED-facet partial trace-mass weights, aligned with the + local ring order. Because every facet is counted exactly once globally, + plain sums of (weight x value) over all ranks' local arrays are exact — + no gather, no ordering, no seam bookkeeping — which is what + :meth:`_surface_mean` reduces over, in any dimension.""" + from underworld3.utilities.boundary_flux import _boundary_stratum_is + + rc = self._ring_coords # local nodes, x-sorted order + if rc.shape[0] == 0 and uw.mpi.size == 1: + return np.empty(0) + dm = self.mesh.dm + cdim = self.mesh.cdim + csec = dm.getCoordinateSection() + cvec = dm.getCoordinatesLocal().array.reshape(-1, cdim) + v0, v1 = dm.getDepthStratum(0) + fS, fE = dm.getHeightStratum(1) + # OWNED facets only: a ghost facet's contribution belongs to the owning + # rank, and each facet must be counted exactly once globally — the seam + # weight is then assembled by SUMMING the per-rank partial contributions + # in the gather (op="sum"), not by averaging copies (a seam copy only + # carries the facets its rank owns). + ghosts = set() + if uw.mpi.size > 1: + nroots, ilocal, _ = dm.getPointSF().getGraph() + if nroots > 0 and ilocal is not None: + ghosts = set(int(i) for i in ilocal) + acc = {} # vertex coord-row -> weight + sis = _boundary_stratum_is(dm, self.mesh, self.surface) + if sis and sis.getSize() > 0: + for f in sis.getIndices(): + if not (fS <= int(f) < fE) or int(f) in ghosts: + continue + verts = [int(p) for p in dm.getTransitiveClosure(int(f))[0] + if v0 <= p < v1] + crows = [csec.getOffset(v) // cdim for v in verts] + xs = cvec[crows] + if len(verts) == 2: # 2D: boundary edge + measure = float(np.linalg.norm(xs[1] - xs[0])) + else: # 3D: boundary triangle + measure = 0.5 * float(np.linalg.norm( + np.cross(xs[1] - xs[0], xs[2] - xs[0]))) + wv = measure / len(verts) + for cr in crows: + acc[cr] = acc.get(cr, 0.0) + wv + # align to the local ring order by position (both sides read the SAME live + # plex coordinates, so the rounded keys match exactly) + wmap = {tuple(np.round(cvec[cr], 9)): w for cr, w in acc.items()} + return np.array([wmap.get(tuple(np.round(c, 9)), 0.0) for c in rc]) def _surface_mean(self, values): r"""Area-weighted global mean of a surface-node array, identical on every rank @@ -368,10 +435,18 @@ def _surface_mean(self, values): The same weighting makes the ``h`` / ``h_inf`` datum volume-preserving rather than node-count-preserving. """ - weights = self._ring_weights() - gathered = self._ring_gather(values) - total = float(weights.sum()) - return float((gathered * weights).sum() / total) if total else 0.0 + # Reduction over OWNED partial weights: each facet's contribution is + # counted exactly once globally, so two scalar allreduces give the exact + # weighted mean with no gather and no ordering — dimension-general (the + # ordered ring remains only for the 2D-only filter and transport). + w = self._surface_weights_local() + v = np.asarray(values, dtype=float) + local_wv = float(np.dot(w, v)) if w.size else 0.0 + local_w = float(w.sum()) if w.size else 0.0 + comm = uw.mpi.comm + total_wv = comm.allreduce(local_wv, op=MPI.SUM) + total_w = comm.allreduce(local_w, op=MPI.SUM) + return total_wv / total_w if total_w else 0.0 def _demean(self, values): """Remove the global surface mean (topography datum floats).""" @@ -520,6 +595,25 @@ def _build_consistent(self): # conds datum, re-evaluated at the boundary nodes at each solve. self.consistent.add_rotated_freeslip_bc( self._un_target.sym[0], self.surface, normal=self.normal) + # Flux-consistent datum gauge: the divergence rows enforce + # ∮ u·n̂_facet over the DEFORMED faceted surface, while the constraint + # fixes u·n̂_node — the directions differ on a deformed surface, so a + # nodal demean of the datum leaves a small net volume flux that an + # incompressible interior cannot absorb (measured: rel ~2e-3 residual + # floor sitting 100% in the pressure rows). Strip the datum's DIRECTED + # mean using the same FE surface integral the residual uses: + # Φ = ∮ ũ_n (n̂·Γ̂) ds and S = ∮ (n̂·Γ̂) ds with Γ̂ = mesh.Gamma (the + # facet normal at quadrature points), then shift ũ_n by Φ/S so the + # discrete net flux of the constrained field is exactly zero. + nrm = (self.normal if self.normal is not None + else self.mesh.boundary_normal(self.surface)) + ncomps = sympy.flatten(sympy.Matrix(nrm)) + ndotg = sum(ncomps[k] * self.mesh.Gamma[k] for k in range(self.mesh.dim)) + self._datum_flux = uw.maths.BdIntegral( + mesh=self.mesh, fn=self._un_target.sym[0] * ndotg, + boundary=self.surface) + self._datum_flux_scale = uw.maths.BdIntegral( + mesh=self.mesh, fn=ndotg, boundary=self.surface) else: n_hat = (self.normal if self.normal is not None else self.mesh.boundary_normal(self.surface)) @@ -561,6 +655,16 @@ def _build_interior_diffuser(self): self._diffuser.constitutive_model = uw.constitutive_models.DiffusionModel self._diffuser.constitutive_model.Parameters.diffusivity = 1.0 self._diffuser.tolerance = 1.0e-3 + # Row map for the deform read: the carrier is P1 and (on our P1-geometry + # meshes) its nodes coincide with the mesh coordinate nodes, so the + # displacement is a DIRECT nodal read — never a point evaluation at the + # field's own nodes (the on-node location class, O(1)-wrong on 3D cell + # edges, #432). Row identities survive deformation; matched once here. + tree = uw.kdtree.KDTree(np.ascontiguousarray(self._carry.coords)) + dist, rows = tree.query(np.ascontiguousarray(self.mesh.X.coords), k=1) + self._carry_rows_at_mesh_nodes = ( + np.asarray(rows).flatten() + if float(np.max(dist)) < 1.0e-12 else None) # exotic geometry: fall back self._base = self._opposite_boundary() self._diffuser.add_essential_bc(self._carry_bc.sym, self.surface) if self._base is not None: @@ -625,19 +729,32 @@ def _build_ring_gather(self): comm = uw.mpi.comm rc = self._ring_coords # local nodes, x-sorted order if self._cylindrical: - s_local = np.arctan2(rc[:, 1], rc[:, 0]) + keys_local = np.arctan2(rc[:, 1], rc[:, 0])[:, None] self._ring_period = 2.0 * np.pi + elif self.mesh.dim == 3: + # No natural 1-D surface ordering in 3D. The ordered-ring FEATURES + # (Taubin filter, tangential transport) are 2D-only and refused at + # construction; the gather itself only needs a deterministic global + # order with exact same-node grouping, which lexicographic rounded + # coordinates provide (the surface deforms along the normal, so the + # reference ordering is built once, like the 2D ring). + keys_local = np.round(np.asarray(rc, dtype=float), 12) + self._ring_period = None else: - s_local = rc[:, 0].astype(float) + keys_local = rc[:, 0].astype(float)[:, None] self._ring_period = None - self._s_local_n = int(s_local.size) + self._s_local_n = int(keys_local.shape[0]) counts = comm.allgather(self._s_local_n) self._ring_offset = int(np.sum(counts[: comm.rank])) - s_global = (np.concatenate(comm.allgather(s_local)) - if comm.size > 1 else s_local.copy()) - self._ring_order = np.argsort(s_global, kind="stable") # concat -> s-sorted + keys_global = (np.concatenate(comm.allgather(np.ascontiguousarray(keys_local))) + if comm.size > 1 else keys_local.copy()) + keys_global = keys_global.reshape(self._s_local_n if comm.size == 1 + else -1, keys_local.shape[1]) + # lexicographic stable order over the key columns (a single column in 2D) + self._ring_order = np.lexsort(keys_global.T[::-1]) self._ring_inv = np.empty_like(self._ring_order) self._ring_inv[self._ring_order] = np.arange(self._ring_order.size) + s_global = keys_global[:, 0] if keys_global.shape[1] == 1 else None # DEDUPLICATE partition-seam copies (#421). A vertex on a partition cut appears # once per adjacent rank in the gathered ring. Every ring operation must see each # PHYSICAL node exactly once: the Taubin filter's roll-stencil otherwise treats @@ -646,31 +763,36 @@ def _build_ring_gather(self): # value — measured as a ~50% seam disagreement in h_inf after 20 iterations and # a few-percent net flux in the prescribed datum. Gather averages the copies # (identical up to round-off); scatter expands back to every copy. - s_srt = s_global[self._ring_order] - uniq_of_sorted = np.empty(s_srt.size, dtype=int) - uniq_id = -1 - prev = None - for i, val in enumerate(np.round(s_srt, 12)): - if prev is None or val != prev: - uniq_id += 1 - prev = val - uniq_of_sorted[i] = uniq_id + keys_srt = np.round(keys_global[self._ring_order], 12) + if keys_srt.shape[0] == 0: + new_group = np.empty(0, dtype=bool) + else: + new_group = np.r_[True, np.any(np.diff(keys_srt, axis=0) != 0.0, axis=1)] + uniq_of_sorted = np.cumsum(new_group) - 1 if keys_srt.shape[0] else \ + np.empty(0, dtype=int) self._ring_uniq_of_sorted = uniq_of_sorted - self._ring_n_uniq = uniq_id + 1 + self._ring_n_uniq = int(uniq_of_sorted[-1] + 1) if keys_srt.shape[0] else 0 self._ring_dup_count = np.bincount(uniq_of_sorted, minlength=self._ring_n_uniq) first_pos = np.searchsorted(uniq_of_sorted, np.arange(self._ring_n_uniq)) - self._s_sorted = s_srt[first_pos] - - def _ring_gather(self, local_vals): - """Local (x-sorted-order) surface values -> the globally s-sorted UNIQUE ring - (partition-seam copies averaged — they agree to round-off by construction).""" + # the 1-D along-surface coordinate exists only where the ordering is real + # (2D); the 3D gather is order-agnostic and the ring features that read + # _s_sorted are refused at construction there. + self._s_sorted = (s_global[self._ring_order][first_pos] + if s_global is not None else None) + + def _ring_gather(self, local_vals, op="mean"): + """Local (x-sorted-order) surface values -> the globally s-sorted UNIQUE ring. + ``op="mean"`` averages partition-seam copies (field values: the copies agree + to round-off by construction); ``op="sum"`` accumulates them (per-rank + PARTIAL contributions such as owned-facet quadrature weights, where the + copies deliberately each carry only their rank's share).""" comm = uw.mpi.comm v = (np.concatenate(comm.allgather(np.ascontiguousarray(local_vals))) if comm.size > 1 else np.asarray(local_vals, dtype=float)) v_sorted = v[self._ring_order] sums = np.bincount(self._ring_uniq_of_sorted, weights=v_sorted, minlength=self._ring_n_uniq) - return sums / self._ring_dup_count + return sums if op == "sum" else sums / self._ring_dup_count def _ring_scatter(self, v_uniq): """Unique-ring values -> this rank's local (x-sorted-order) nodes. Every seam @@ -960,6 +1082,16 @@ def _solve_consistent(self, increment, dt): u_tilde = self._demean(increment / dt) self._un_target.array[...] = 0.0 self._un_target.array[self._un_target_rows, 0, 0] = u_tilde + # Exact FE-consistent flux strip (STRONG constraint only — the penalty + # absorbs a datum flux weakly and builds no integrals): remove the + # DIRECTED mean so the datum carries zero discrete net flux through the + # deformed facets — the quantity the pressure rows actually enforce. + # Collective (BdIntegral), so the shift is identical on every rank. + if self.consistent_constraint == "strong": + flux = float(self._datum_flux.evaluate()) + scale = float(self._datum_flux_scale.evaluate()) + if abs(scale) > 1.0e-30: + self._un_target.array[self._un_target_rows, 0, 0] -= flux / scale # Warm-start from the free solve: the consistent solution IS the free # solution with the (small) material-boundary datum imposed, and the free # solve has already converged this step. Starting there keeps a power-law @@ -990,8 +1122,15 @@ def _carry_and_deform(self, increment, dt): self._carry_bc.array[self._carry_bc_rows, 0, 0] = increment self._diffuser.solve(zero_init_guess=False) coords = self.mesh.X.coords - displacement = np.asarray( - function.evaluate(self._carry.sym[0], coords) - ).flatten() + if self._carry_rows_at_mesh_nodes is not None: + # direct nodal read (see _build_interior_diffuser: the on-node + # evaluation class is what this avoids) + displacement = np.asarray( + self._carry.array[self._carry_rows_at_mesh_nodes, 0, 0] + ).flatten() + else: + displacement = np.asarray( + function.evaluate(self._carry.sym[0], coords) + ).flatten() new_coords = coords + displacement[:, None] * self._normal_direction(coords) self.mesh.deform(new_coords, dt=dt) diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index 4379c934..6142ec69 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/src/underworld3/utilities/boundary_flux.py @@ -165,8 +165,9 @@ def _desmear(solver, boundary, xs, R, mass, remove_mean, partial_reaction=True): csec = dm.getCoordinateSection() cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) v0, v1 = dm.getDepthStratum(0) - if mass not in ("auto", "lumped", "consistent"): - raise ValueError("mass must be 'auto', 'lumped', or 'consistent'.") + if mass not in ("auto", "lumped", "consistent", "p1"): + raise ValueError("mass must be 'auto', 'lumped', 'consistent', or 'p1' " + "(P1-projected recovery on a 3D P2 trace).") if dim == 3: lsec = dm.getLocalSection() ncomp = lsec.getFieldComponents(0) @@ -260,8 +261,44 @@ def coord(q): if order == 2 and mass == "lumped": raise ValueError( "A 3D P2 triangular trace has zero row-sum mass at its vertices; " - "use mass='consistent' for pointwise boundary-flux recovery." + "use mass='consistent' (pointwise, carries the vertex-integral " + "checkerboard risk) or mass='p1' (P1-projected, monotone — the " + "choice for driving a P1 surface field) for boundary-flux recovery." ) + mid_owners = {} + if mass == "p1": + if order != 2: + mass = "lumped" # P1 trace: p1 IS lumped + else: + # P1-PROJECTED recovery on a P2 trace: the consistent P2 path has + # the ∫φ_vertex = 0 vertex checkerboard (#404 hold), while the P1 + # trace is sound — and a P1 surface field only consumes vertex + # values anyway. Fold each edge-midpoint load onto its two edge + # vertices (φ^{P1}(edge-mid) = 1/2 exactly, P1 ⊂ P2 — the load + # transfer is the interpolation transpose, so the total load is + # conserved), then de-smear with the P1 lumped triangle mass. + # Midpoint outputs are read back as the P1 interpolant (vertex + # average). + new_elements = {} + for _order, nodes, area in elements.values(): + vk = nodes[:3] + m01, m12, m20 = nodes[3:] + mid_owners[m01] = (vk[0], vk[1]) + mid_owners[m12] = (vk[1], vk[2]) + mid_owners[m20] = (vk[2], vk[0]) + new_elements[(1, tuple(sorted(vk)))] = (1, vk, area) + folded = {} + for key, value in R_by.items(): + if key in mid_owners: + va, vb = mid_owners[key] + folded[va] = folded.get(va, 0.0) + 0.5 * value + folded[vb] = folded.get(vb, 0.0) + 0.5 * value + else: + folded[key] = folded.get(key, 0.0) + value + R_by = folded + elements = new_elements + order = 1 + mass = "lumped" keys = sorted(R_by) global_index = {key: i for i, key in enumerate(keys)} @@ -311,14 +348,23 @@ def coord(q): if remove_mean: mean = float(np.dot(flux, boundary_mass) / np.sum(boundary_mass)) flux -= mean - return np.array([flux[global_index[_key(x, dim)]] for x in xs]) + + def value_at(x): + key = _key(x, dim) + if key in global_index: + return flux[global_index[key]] + # P1-projected mode: a P2 edge midpoint reads the P1 interpolant + va, vb = mid_owners[key] + return 0.5 * (flux[global_index[va]] + flux[global_index[vb]]) + + return np.array([value_at(x) for x in xs]) if dim != 2: raise NotImplementedError( f"Boundary-flux recovery is not implemented for mesh dimension {dim}." ) - if mass == "auto": - mass = "lumped" + if mass in ("auto", "p1"): + mass = "lumped" # 2D: the lumped line mass is sound e0, e1 = dm.getDepthStratum(1) def vcoord(q): return cvec[csec.getOffset(q) // dim] diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 448a7297..19f61f80 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -645,6 +645,19 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, diag_scale = None lin_its = [] + # With the constant-pressure nullspace active, the outer residual is measured + # in the pressure-gauge QUOTIENT space: the inner KSP projects the constant + # mode out of every increment (it is the gauge of an enclosed incompressible + # domain), so the loop cannot reduce that component and must not measure it. + # On a DEFORMED faceted boundary the component is not exactly zero — free + # tangential DOFs carry a small net flux through the node-vs-facet normal + # mismatch, an irreducible discrete incompatibility (measured: an unprojected + # outer norm floors at rel ~2e-3, 100% pressure rows, and the line search + # stalls against the constant offset). This mirrors PETSc's own projected + # residual for singular systems. The Cartesian reaction stash is UNPROJECTED + # (σ_nn reads velocity rows only). + use_pnull = bool(getattr(solver, "_petsc_use_pressure_nullspace", False)) + def rotated_residual(uvec, keep_cartesian=False): snes.computeFunction(uvec, Fc) if keep_cartesian: @@ -652,6 +665,11 @@ def rotated_residual(uvec, keep_cartesian=False): Fh = Fc.duplicate() Q.mult(Fc, Fh) _zero_rows_local(Fh, normal_rows) + if use_pnull: + sp = Fh.getSubVector(pres_is) + mean = sp.sum() / max(sp.getSize(), 1) + sp.shift(-mean) # project out the pressure-gauge mode + Fh.restoreSubVector(pres_is, sp) return Fh # Convergence reference: max(initial residual, REST-STATE residual ‖F̂(0)‖). diff --git a/tests/test_1070_free_surface_plume.py b/tests/test_1070_free_surface_plume.py index 8904b70b..a0d6b441 100644 --- a/tests/test_1070_free_surface_plume.py +++ b/tests/test_1070_free_surface_plume.py @@ -311,8 +311,19 @@ def test_freesurface_ring_quadrature_is_exact_in_parallel(): fs = uw.systems.FreeSurface(stokes, "Upper", buoyancy_scale=1.0, normal=rhat) weights = fs._ring_weights() - assert abs(float(weights.sum()) - 2.0 * np.pi * r_out) < 1.0e-9, \ - f"ring weights sum to {weights.sum():.10f}, not the circumference (seam double-count?)" + # The gauge must match the FE boundary quadrature EXACTLY: the datum's + # flux-free condition lives in the discrete space, so the weight total is the + # FE measure of the (polygonal) boundary — NOT the ideal-circle arc length, + # which differs at O(h^2) and was the strong-datum compatibility floor + # (block-split evidence: the whole stalled residual sat in the pressure rows). + # A seam double-count breaks this equality at the first shared node. + fe_len = float(uw.maths.BdIntegral(mesh=mesh, fn=sympy.S.One, + boundary="Upper").evaluate()) + assert abs(float(weights.sum()) - fe_len) < 1.0e-9, \ + f"ring weights sum to {weights.sum():.10f}, FE boundary measure is {fe_len:.10f} " \ + "(seam double-count?)" + # sanity: the polygonal measure approximates the circle at O(h^2) + assert abs(fe_len - 2.0 * np.pi * r_out) < 1.0e-2 coords = fs._ring_coords theta = np.arctan2(coords[:, 1], coords[:, 0]) diff --git a/tests/test_1072_free_surface_spherical.py b/tests/test_1072_free_surface_spherical.py new file mode 100644 index 00000000..7e83d48a --- /dev/null +++ b/tests/test_1072_free_surface_spherical.py @@ -0,0 +1,59 @@ +"""3D FreeSurface on a spherical shell: the end-to-end loop must run and produce +physically sensible topography (guards the dimension-general surface machinery: +owned-facet trace-mass gauge, P1-projected sigma_nn recovery, radial deform). + +The quantitative benchmarking (analytic Y_lm shell rate, convergence of the +h_inf modal bias, 3D parallel) is the review-team's scope — this test pins the +CAPABILITY: construction, one solve/advance cycle, finite mean-free h_inf, and +the explicit refusal of the 2D-only features. +""" +import numpy as np +import pytest +import sympy +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + + +def _shell_stokes(cell=0.35): + mesh = uw.meshing.SphericalShell(radiusOuter=1.0, radiusInner=0.547, + cellSize=cell, qdegree=3) + x, y, z = mesh.X + r = sympy.sqrt(x ** 2 + y ** 2 + z ** 2) + rhat = sympy.Matrix([[x / r, y / r, z / r]]) + stokes = uw.systems.Stokes(mesh) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + blob = sympy.exp(-(((x - 0.75) ** 2 + y ** 2 + z ** 2) / 0.05)) + stokes.bodyforce = 50.0 * blob * rhat.T + stokes.add_essential_bc((0.0, 0.0, 0.0), "Lower") + stokes.tolerance = 1.0e-5 + return mesh, stokes, rhat + + +def test_freesurface_spherical_shell_end_to_end(): + """Construction + one full solve/advance on the shell; h_inf finite and + mean-free; the surface responds toward equilibrium (|h| grows from flat + under the one-sided load and stays bounded by |h_inf|).""" + mesh, stokes, rhat = _shell_stokes() + fs = uw.systems.FreeSurface(stokes, "Upper", buoyancy_scale=50.0, normal=rhat) + fs.solve() + h_inf = np.asarray(fs._h_inf) + assert np.isfinite(h_inf).all(), "3D h_inf recovery produced non-finite values" + assert abs(fs._surface_mean(h_inf)) < 1.0e-8 * (np.abs(h_inf).max() + 1e-30), \ + "h_inf datum is not mean-free under the trace-mass gauge" + assert np.abs(h_inf).max() > 1.0e-4, "no topographic response to the load" + fs.advance(fs.estimate_dt(advect_scale=10.0)) + shape = fs._current_shape() + assert np.isfinite(shape).all() + assert 0.0 < np.abs(shape).max() <= 1.5 * np.abs(h_inf).max(), \ + "surface did not move toward (or overshot) equilibrium" + + +def test_freesurface_spherical_refuses_2d_only_features(): + """The 2D-only features fail loudly at construction in 3D, not silently.""" + mesh, stokes, rhat = _shell_stokes() + with pytest.raises(NotImplementedError, match="tangential"): + uw.systems.FreeSurface(stokes, "Upper", normal=rhat, tangent_advect="shape") + with pytest.raises(NotImplementedError, match="filter"): + uw.systems.FreeSurface(stokes, "Upper", normal=rhat, surface_filter=10)