diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index f56fbe0fe..8ab13743e 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -4315,7 +4315,7 @@ class SNES_Vector(SolverBaseClass): def add_nitsche_bc(self, conds=None, boundary=None, direction=None, - normal=None, gamma=10.0, theta=1, mask=None, + normal=None, gamma=12.5, theta=1, mask=None, local_h=True, g=None): r"""Add Nitsche weak enforcement of a velocity constraint along a direction. @@ -4337,8 +4337,13 @@ class SNES_Vector(SolverBaseClass): terms — the same geometric-normal override as on the Stokes variant. Default ``None`` uses the per-boundary, deformation-tracking ``mesh.boundary_normal(boundary)``. - gamma : float, default=10.0 - Dimensionless stabilisation parameter. + gamma : float, default=12.5 + Dimensionless stabilisation parameter. The penalty is + ``gamma*mu/h``, so this is calibrated against the definition of + ``h``. It was 10.0 while ``h`` came from a kd-tree of neighbouring + centroids; ``mesh.cell_size()`` is now PETSc's ``volume**(1/dim)``, + about 19% larger in the mean, and 12.5 restores the enforcement + that gamma=10 gave against the old h (#694). theta : {-1, 0, 1}, default=1 Symmetry parameter (1=symmetric, -1=skew-symmetric). mask : sympy expression, optional @@ -6532,7 +6537,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): remove_mean=remove_mean) def add_nitsche_bc(self, conds=None, boundary=None, direction=None, normal=None, - gamma=10.0, theta=1, mask=None, local_h=True, g=None): + gamma=12.5, theta=1, mask=None, local_h=True, g=None): r"""Add Nitsche weak enforcement of a velocity constraint along a direction. Nitsche's method provides a variationally consistent alternative to @@ -6569,9 +6574,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): Boundary unit normal used in the Nitsche consistency, symmetry, and pressure-coupling terms. Default ``None`` uses the per-boundary, deformation-tracking ``mesh.boundary_normal(boundary)``. - gamma : float, default=10.0 + gamma : float, default=12.5 Dimensionless stabilisation parameter. Typical values 5--20 - for P2 elements. + for P2 elements. The penalty is ``gamma*mu/h``, so this is + calibrated against the definition of ``h``: it was 10.0 while + ``h`` came from a kd-tree of neighbouring centroids, and moved + with ``mesh.cell_size()`` becoming PETSc's ``volume**(1/dim)`` + (#694). theta : {-1, 0, 1}, default=1 Symmetry parameter: 1: symmetric (default — optimal convergence and solver efficiency) diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 56a5222a0..e3963824a 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/src/underworld3/discretisation/discretisation_mesh.py @@ -2840,12 +2840,10 @@ def nuke_coords_and_rebuild( flush=True, ) - ( - self._min_size, - self._radii, - self._centroids, - self._search_lengths, - ) = self._get_mesh_sizes() + # `_min_size` and `_search_lengths` used to be unpacked here and were + # never read anywhere in src/ or tests/ -- the kd-tree loop computed + # three distance statistics per cell and two were discarded. + self._cell_radii, self._centroids = self._get_cell_radii() # Skip self-copy when hierarchy is trivial (issue #96 investigation) if self.dm is not self.dm_hierarchy[-1]: @@ -3234,8 +3232,8 @@ def cell_size(self): Returns the ``.sym`` of a cell-constant (degree-0, discontinuous) scalar MeshVariable holding each cell's characteristic length (the - RMS distance of its vertices from their own centroid). This is a - purely cell-local quantity, independent of the MPI partition. Unlike + ``volume**(1/dim)`` equivalent radius, i.e. ``self._cell_radii``, + which comes from PETSc and is independent of the MPI partition). Unlike the single *global* scalar from :meth:`get_min_radius` (the smallest cell anywhere), this varies cell to cell, so a stabilisation that scales as :math:`1/h` — e.g. the Nitsche free-slip penalty @@ -3301,27 +3299,45 @@ def _refresh(): def _assemble_cell_size(self, var): """Fill ``var`` (degree-0 scalar) with each cell's characteristic size. - Uses the cell-geometry characteristic lengths ``self._cell_radii`` computed by - :meth:`_get_mesh_sizes` on the *current* geometry. A degree-0 + Uses the per-cell characteristic lengths ``self._cell_radii`` computed + by :meth:`_get_cell_radii` on the *current* geometry. A degree-0 discontinuous variable's local DOFs and ``self._cell_radii`` are BOTH indexed by this rank's cell-stratum order, so a direct assignment is correct on every rank. - This is deliberately a purely RANK-LOCAL operation (no ``var.coords`` - access, no collective): mixing a rank-local fast path with a - collective fallback would diverge across ranks and deadlock, because - ``var.coords`` triggers the collective ``_get_coords_for_basis``.""" - # Own-cell radii fix #687 without changing the legacy kd-tree radii - # used by global timestep estimates, adaptivity, and mesh relaxation. + This routine has TWO collectives in it and therefore no early return, + which is the opposite of what its previous docstring claimed (#698): + + * the first ``var.data`` access lazily reaches ``MeshVariable._set_vec``, + which calls ``dm.createSubDM`` and ``createGlobalVector``; + * assigning into ``var.data`` fires the array's write-back callback, + ``pack_raw_data_to_petsc``. + + A rank owning no cells used to short-circuit on ``radii.size == 0`` and + return before either, while its populated peers made both. Measured on a + region submesh at np=8 with cells per rank ``[12, 11, 0, 19, 0, 0, 0, 0]``: + the populated ranks sat in ``pack_raw_data_to_petsc`` and the job never + finished. A starved rank now walks the same path writing a zero-length + slice. + + ``var.coords`` is still deliberately not read: it triggers the collective + ``_get_coords_for_basis``, and reading it only on some ranks would put a + third conditional collective back in.""" + # `_cell_radii` is PETSc's volume**(1/dim), a property of each cell, so + # the values here do not depend on the partition -- and neither does the + # Nitsche penalty gamma*mu/h that consumes them under the default + # local_h=True. It was a kd-tree distance to the nearest centroid among + # THIS RANK's centroids, which near a seam could simply be absent (#694). + # There is NO early return here, and that is the point. Both steps below + # are collective, so a rank owning no cells has to walk through them + # writing nothing rather than skipping them (#698). + data = var.data # allocates: createSubDM + createGlobalVector radii = numpy.asarray(self._cell_radii).reshape(-1) - # Empty partition (no local cells): nothing to fill on this rank. - if radii.size == 0 or var.data.shape[0] == 0: - return - # Assign over the common length. In practice these match exactly (same - # local cell set / ordering); the slice only guards a stray off-by-ghost - # mismatch without ever taking a collective path on a subset of ranks. - n = min(var.data.shape[0], radii.shape[0]) - var.data[:n, 0] = radii[:n] + + # `n` is 0 on a starved rank. The assignment still fires the array's + # write-back callback, which is the second collective. + n = min(data.shape[0], radii.shape[0]) + data[:n, 0] = radii[:n] @property def Gamma_P1(self): @@ -6889,56 +6905,44 @@ def _eval_use_robust_location(self) -> bool: """ return (uw.mpi.size > 1) and (self._location_capability() != "none") - def _get_mesh_sizes(self, verbose=False): - """ - Cache own-cell radii for cell_size and return legacy kd-tree radii. - - Own-cell sizes use current DM vertices, so neither partition-local - neighbours nor stale coordinate views affect stabilization (#687). - Legacy radii remain unchanged for their other consumers. - """ - + def _get_cell_radii(self): + """Each cell's characteristic length, and the cell centroids. + + The length is PETSc's ``volume**(1/dim)`` from + ``DMPlexComputeGeometryFVM``. A cell's volume is a property of that + cell, so this cannot depend on how the mesh was partitioned — which is + the point. + + It replaces a kd-tree of THIS RANK's centroids queried with each cell's + vertices. Near a partition boundary the true nearest centroid can belong + to a cell owned by another rank and be absent from the tree, so the + answer moved with the rank count: per-cell by 3.3e-03 at np=2 and + 4.1e-03 at np=4, `get_max_radius()` by 4.9% at np=8, `get_mean_radius()` + at every rank count, and `mesh.cell_size()` with them -- which scales + the Nitsche penalty under the DEFAULT ``local_h=True`` (#569, #687, + #694). + + The FVM routine had been abandoned with a note that it "does not + compute all cells". That does not reproduce: measured on 2-D simplex, + 2-D quad, 3-D tetrahedra, 3-D hexahedra and a deformed mesh, it returns + one finite positive value per local cell and is bit-identical across + rank counts in every case. (The note also named ``DMPlexGetMinRadius``, + which is a different call and is not used here.) + """ + from underworld3.cython import petsc_discretisation + + radii, _fvm_centroids = petsc_discretisation.petsc_fvm_get_local_cell_sizes(self) + + # The FVM centroids are discarded: `_get_coords_for_basis(0, False)` is + # the degree-0 coordinate array the rest of the mesh indexes by cell, + # and mixing the two orderings would misalign every per-cell lookup. centroids = self._get_coords_for_basis(0, False) - centroids_kd_tree = uw.kdtree.KDTree(centroids) - import numpy as np - - cStart, cEnd = self.dm.getHeightStratum(0) - pStart, pEnd = self.dm.getDepthStratum(0) - cell_length = np.empty(centroids.shape[0]) - cell_min_r = np.empty(centroids.shape[0]) - cell_r = np.empty(centroids.shape[0]) - cell_radii = np.empty(centroids.shape[0]) - coordinate_section = self.dm.getCoordinateDM().getLocalSection() - vertex_coordinates = self.dm.getCoordinatesLocal().array - - for cell in range(cEnd - cStart): - cell_num_points = self.dm.getConeSize(cell) - cell_points = self.dm.getTransitiveClosure(cell)[0][-cell_num_points:] - # Use raw internal array for internal mesh operations (avoid unit-aware wrapping) - cell_coords = self._coords[cell_points - pStart] - - distsq, _ = centroids_kd_tree.query(cell_coords, k=1, sqr_dists=True) - - cell_length[cell] = np.sqrt(distsq.max()) - cell_r[cell] = np.sqrt(distsq.mean()) - cell_min_r[cell] = np.sqrt(distsq.min()) - - # A hex has six faces but eight vertices: select the vertex - # stratum, not a cone-sized suffix of its transitive closure. - closure = self.dm.getTransitiveClosure(cStart + cell)[0] - vertices = closure[(closure >= pStart) & (closure < pEnd)] - offsets = np.array([coordinate_section.getOffset(int(v)) for v in vertices]) - own_coords = vertex_coordinates[offsets[:, None] + np.arange(self.cdim)] - delta = own_coords - own_coords.mean(axis=0) - cell_radii[cell] = np.sqrt(np.mean(np.sum(delta ** 2, axis=1))) - - self._cell_radii = cell_radii - return cell_min_r, cell_r, centroids, cell_length + return radii, centroids # ========== - # Deprecated in favour of _get_mesh_sizes (above) + # Deprecated in favour of _get_cell_radii (above) def _get_mesh_centroids(self): """ Obtain and cache the (local) mesh centroids using underworld swarm technology. @@ -7028,7 +7032,7 @@ def get_min_radius(self) -> float: import numpy as np from mpi4py import MPI - radii = np.asarray(self._radii).reshape(-1) + radii = np.asarray(self._cell_radii).reshape(-1) local_min = float(radii.min()) if radii.size else float("inf") if uw.mpi.size > 1: local_min = uw.mpi.comm.allreduce(local_min, op=MPI.MIN) @@ -7050,7 +7054,7 @@ def get_max_radius(self) -> float: import numpy as np from mpi4py import MPI - radii = np.asarray(self._radii).reshape(-1) + radii = np.asarray(self._cell_radii).reshape(-1) local_max = float(radii.max()) if radii.size else float("-inf") if uw.mpi.size > 1: local_max = uw.mpi.comm.allreduce(local_max, op=MPI.MAX) @@ -7069,15 +7073,22 @@ def get_mean_radius(self) -> float: this is the canonical "mesh length" API. Use this anywhere you need a representative h0 (smoothing-length defaults, diffusion- stability heuristics, problem-scale normalisation) rather than - reaching for the rank-local ``self._radii`` array, which gives - different answers on different MPI ranks and leaks downstream - (e.g. into JIT C source via per-rank pointwise-function inputs). + reducing a per-rank array by hand, which gives different answers on + different MPI ranks and leaks downstream (e.g. into JIT C source via + per-rank pointwise-function inputs). + + The value is the same at every RANK COUNT as well as on every rank: + ``self._cell_radii`` is PETSc's ``volume**(1/dim)``, a property of each + cell rather than of the partition. That was not true while these + reduced over a kd-tree of this rank's centroids -- an allreduce made + the answer agree across ranks without making it agree across rank + counts, and ``get_max_radius()`` moved 4.9% at np=8 (#694). """ import numpy as np from mpi4py import MPI - radii = np.asarray(self._radii) + radii = np.asarray(self._cell_radii) local_sum = float(radii.sum()) local_n = int(radii.size) if uw.mpi.size > 1: diff --git a/src/underworld3/meshing/smoothing/api.py b/src/underworld3/meshing/smoothing/api.py index 145e7977d..5b4160085 100644 --- a/src/underworld3/meshing/smoothing/api.py +++ b/src/underworld3/meshing/smoothing/api.py @@ -713,7 +713,7 @@ def follow_metric( mesh, T, refinement=2.0, coarsening=2.0, metric="gradient-uniform", - gradient_smoothing_length=2.0 * mesh._radii.mean(), + gradient_smoothing_length=2.0 * mesh.get_mean_radius(), ) See Also diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index 4b4b679d3..45ca3f18b 100644 --- a/src/underworld3/swarm.py +++ b/src/underworld3/swarm.py @@ -5075,7 +5075,7 @@ def estimate_dt(self, V_fn): vel = uw.function.evaluate(V_fn, self._particle_coordinates.data, evalf=True) # If vel is unit-aware (UnitAwareArray), nondimensionalise it to get - # consistent nondimensional values that match mesh._radii + # consistent nondimensional values that match mesh._cell_radii # Note: .magnitude returns physical units, which would be wrong here if hasattr(vel, "units") and vel.units is not None: vel = uw.non_dimensionalise(vel) diff --git a/src/underworld3/systems/free_surface.py b/src/underworld3/systems/free_surface.py index 6fd131772..dc6f803c9 100644 --- a/src/underworld3/systems/free_surface.py +++ b/src/underworld3/systems/free_surface.py @@ -121,7 +121,7 @@ def __init__( driving_buoyancy=None, smooth_length=0.0, mass="lumped", - max_surface_cfl=0.5, + max_surface_cfl=0.437, tangent_advect=None, tangent_spectral_modes=0, surface_mask=None, diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index c89963f99..54fe7a120 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -256,7 +256,7 @@ def _global_max_diffusivity(constitutive_K, mesh): diffusivity = K # If unit-aware (UnitAwareArray), nondimensionalise so the value is - # consistent with mesh._radii. Note: .magnitude alone would keep the + # consistent with mesh._cell_radii. Note: .magnitude alone would keep the # physical-units number, which would be wrong here. if hasattr(diffusivity, "units") and diffusivity.units is not None: diffusivity = uw.non_dimensionalise(diffusivity) @@ -276,7 +276,7 @@ def _centroid_velocities_nd(V_fn, mesh, basis=None, ensure_2d=True): Shared by the ``estimate_dt`` implementations: the advective CFL limit needs per-element centroid velocities in the same (nondimensional) - scale as ``mesh._radii``. + scale as ``mesh._cell_radii``. Parameters ---------- @@ -303,7 +303,7 @@ def _centroid_velocities_nd(V_fn, mesh, basis=None, ensure_2d=True): vel = uw.function.evaluate(V_fn, mesh._centroids) # If unit-aware (UnitAwareArray), nondimensionalise so the values are - # consistent with mesh._radii. Note: .magnitude alone would keep the + # consistent with mesh._cell_radii. Note: .magnitude alone would keep the # physical-units numbers, which would be wrong here. if hasattr(vel, "units") and vel.units is not None: vel = uw.non_dimensionalise(vel) @@ -2319,7 +2319,7 @@ def estimate_dt(self): vel_magnitudes = np.linalg.norm(vel, axis=1) # Get per-element radii (characteristic element size) - element_radii = self.mesh._radii + element_radii = self.mesh._cell_radii # Compute per-element advective timestep: dt_i = h_i / |v_i| # Avoid division by zero for elements with zero velocity @@ -4336,7 +4336,7 @@ def estimate_dt(self, direction_aware: bool = False, percentile: float = 0.0): centroid) · v̂` over the cell vertices. This is the distance material actually traverses through the cell per unit ``|v|``, and is **always ≥ the isotropic - mesh._radii estimate**, by 1.5–3× for equant cells + mesh._cell_radii estimate**, by 1.5–3× for equant cells (geometric factor) and up to ~10× for cells that the mover has stretched along the flow direction. On adapted meshes the gain is substantial; on uniform @@ -4379,7 +4379,7 @@ def estimate_dt(self, direction_aware: bool = False, percentile: float = 0.0): vel_magnitudes = np.linalg.norm(vel, axis=1) # Get per-element radii (characteristic element size) - element_radii = self.mesh._radii + element_radii = self.mesh._cell_radii ## estimate dt of adv and diff components using per-element approach ## dt_adv_i = h_i / |v_i| for advection @@ -4409,7 +4409,7 @@ def _reduce_dt(per_elem): dt_diff_per_element = np.array([np.inf]) # Per-element advective timestep — either isotropic - # (mesh._radii / |v|) or direction-aware (v-aligned cell + # (mesh._cell_radii / |v|) or direction-aware (v-aligned cell # extent / |v|). if direction_aware: # Per-cell vertex indices (triangle / tet). diff --git a/tests/parallel/test_0774_empty_rank_reductions_mpi.py b/tests/parallel/test_0774_empty_rank_reductions_mpi.py index e763ddb8d..da93112fe 100644 --- a/tests/parallel/test_0774_empty_rank_reductions_mpi.py +++ b/tests/parallel/test_0774_empty_rank_reductions_mpi.py @@ -1,7 +1,7 @@ """Parallel regression tests for issue #405 — reductions on a zero-cell rank. A rank that owns NO CELLS used to raise a rank-local ``ValueError`` from an -unguarded local reduction (``self._radii.min()`` and friends) while its +unguarded local reduction (``self._cell_radii.min()`` and friends) while its populated peers sat in the matching collective. The job then hung or aborted asymmetrically. Every global quantity computed from rank-local data must instead reduce across ranks, with the starved rank contributing the identity @@ -89,7 +89,7 @@ def test_negative_control_rank_local_minimum_would_be_caught(): """ mesh = _starved_box() - radii = np.asarray(mesh._radii).reshape(-1) + radii = np.asarray(mesh._cell_radii).reshape(-1) rank_local_min = float(radii.min()) if radii.size else float("inf") gathered = uw.mpi.comm.allgather(rank_local_min) diff --git a/tests/parallel/test_0779_submesh_cell_size_starved_rank.py b/tests/parallel/test_0779_submesh_cell_size_starved_rank.py new file mode 100644 index 000000000..80358dcb5 --- /dev/null +++ b/tests/parallel/test_0779_submesh_cell_size_starved_rank.py @@ -0,0 +1,69 @@ +"""``cell_size()`` must survive a rank that owns none of the submesh. + +A region submesh keeps only the cells carrying its label, so a partition whose +share of the parent lies outside that region legitimately owns nothing. That is +routine for submeshes and is the case this guards -- not an over-decomposed +full mesh, which is a different (and degenerate) situation. + +The defect (#698): ``_assemble_cell_size`` returned early on ``radii.size == 0``, +in front of TWO collectives --- the first ``var.data`` access, which allocates +through ``dm.createSubDM`` and ``createGlobalVector``, and the assignment into +it, which fires ``pack_raw_data_to_petsc``. Starved ranks skipped both while +their populated peers made them, and the job never finished. Measured at np=8 +with cells per rank ``[12, 11, 0, 19, 0, 0, 0, 0]``. +""" + +import pytest + +import underworld3 as uw + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_1, pytest.mark.tier_a] + + +def _thin_slab_parent(): + """A box split near the top, so the outer region is a few cells deep.""" + return uw.meshing.BoxInternalBoundary( + elementRes=(12, 12), zelementRes=(10, 2), zintCoord=0.85, + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), simplex=True, qdegree=2) + + +@pytest.mark.mpi(min_size=2) +def test_cell_size_completes_when_ranks_own_none_of_the_submesh(): + """The whole test is that this returns. A regression hangs rather than fails. + + The starved-rank count is reported, not asserted: how the parent is split is + the partitioner's business, and at np=2 the slab may well reach both ranks. + What must hold at every rank count is that every rank comes back. + """ + submesh = _thin_slab_parent().extract_region("Outer") + start, end = submesh.dm.getHeightStratum(0) + counts = uw.mpi.comm.allgather(end - start) + + submesh.cell_size() + uw.mpi.comm.barrier() + + uw.mpi.pprint(f"SUBMESH_CELL_SIZE ranks={uw.mpi.size} cells={counts} " + f"starved={counts.count(0)}") + assert sum(counts) > 0, "the region submesh is empty everywhere; test is vacuous" + + +@pytest.mark.mpi(min_size=2) +def test_the_size_field_is_filled_where_there_are_cells(): + """A starved rank writing nothing must not stop the others writing. + + Without this, the fix above could be satisfied by never filling the field + at all. + """ + import numpy as np + + submesh = _thin_slab_parent().extract_region("Outer") + start, end = submesh.dm.getHeightStratum(0) + submesh.cell_size() + + local = np.asarray(submesh._cell_size_variable.array[:, 0, 0]) + assert local.shape[0] == end - start + if local.size: + assert np.all(local > 0.0), "populated rank has non-positive cell sizes" + from mpi4py import MPI + filled = uw.mpi.comm.allreduce(int(local.size), op=MPI.SUM) + assert filled > 0, "no rank filled the field at all" diff --git a/tests/parallel/test_1069_boundary_normal_parallel.py b/tests/parallel/test_1069_boundary_normal_parallel.py index 036f39a45..152b25101 100644 --- a/tests/parallel/test_1069_boundary_normal_parallel.py +++ b/tests/parallel/test_1069_boundary_normal_parallel.py @@ -333,12 +333,14 @@ def _nitsche_annulus_diagnostics(): leakage. Both are stable from tolerance 1e-9 to 1e-12, so neither is the linear solve. - This test now leaves ``local_h`` at its default ``True``. Before #569/#687, - doing so mixed the boundary-normal regression with a second partition-dependent - input from ``mesh.cell_size()``; this test therefore had to disable the public - default. The cell-local geometric size is now partition independent, so retaining - the default jointly guards the normal assembly and the Nitsche penalty path users - actually run. + This runs with the DEFAULT ``local_h=True``. It used to pass ``local_h=False`` + because ``mesh.cell_size()`` was itself partition-dependent -- a kd-tree over + THIS RANK's cell centroids, so on this mesh the field's sum was 26.0822 at + np=1, 26.1211 at np=2 and 26.1386 at np=4 -- and leaving the default on would + have made this test measure two defects at once. ``cell_size()`` is now + PETSc's ``volume**(1/dim)`` and is partition independent (#569, #687, #694), + so keeping the default guards the normal assembly and the Nitsche penalty + path users actually run, together. """ RI, RO = 0.5, 1.0 mesh = uw.meshing.Annulus(radiusInner=RI, radiusOuter=RO, cellSize=0.12, qdegree=3) diff --git a/tests/parallel/test_1078_radius_accessors_partition_independence.py b/tests/parallel/test_1078_radius_accessors_partition_independence.py new file mode 100644 index 000000000..f3d5a0433 --- /dev/null +++ b/tests/parallel/test_1078_radius_accessors_partition_independence.py @@ -0,0 +1,119 @@ +"""The radius ACCESSORS must not depend on how the mesh was partitioned. + +``test_1077`` covers the per-cell field cell by cell. This file covers +``get_min_radius()``, ``get_max_radius()`` and ``get_mean_radius()``, which +reduce that field and advertise a global mesh length -- and which were the half +left partition-dependent when only ``cell_size()`` was fixed: ``get_max_radius()`` +moved 4.9% at np=8 and ``get_mean_radius()`` at every rank count. + +The defect (#569, #687, #694): the per-cell length came from a kd-tree over +THIS RANK's centroids, queried with each cell's vertices. Near a partition +boundary the true nearest centroid can belong to a cell owned by another rank +and so be absent from the tree, and the answer moved with the rank count -- +per-cell by 3.3e-03 at np=2, ``get_max_radius()`` by 4.9% at np=8, +``get_mean_radius()`` at every rank count. + +It reached users through ``mesh.cell_size()``, which scales the Nitsche penalty +under the default ``local_h=True``, and through the three radius accessors, +whose docstrings advertise a global mesh length. + +Two things this file is careful about, both learned the hard way: + +* **It compares two rank counts.** Partition independence is a statement about + two runs agreeing, so nothing measured at a single rank count establishes it. + The tests that shipped with the first attempt at this fix asserted a + within-rank oracle -- each cell matching its own vertices -- which is true of + a partition-dependent field as well. +* **The reference is computed here, not recorded.** ``serial_reference`` runs + this module's own ``__main__`` at np=1 in this environment and asserts the + mesh fingerprints match, so a host that triangulates differently is reported + as that rather than as partition dependence. +""" + +import numpy as np +import pytest + +import underworld3 as uw +from mpi4py import MPI + +from serial_reference import compare, emit, mesh_fingerprint, serial_reference + +pytestmark = [pytest.mark.mpi(min_size=2), pytest.mark.level_1, pytest.mark.tier_a] + +LABELS = ("min radius", "max radius", "mean radius", "sum of cell radii") + +# min and max are exact reductions of identical per-cell values, so they must +# agree to the bit. The mean is a distributed sum and reduces in partition +# order, so it is allowed the last couple of bits -- that ordering difference +# is not the defect under test. +RTOLS = (0.0, 0.0, 1.0e-12, 1.0e-12) + + +def _box(): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.12, qdegree=2) + # mesh_fingerprint integrates over the mesh, and integration needs at + # least one variable to exist on it. + uw.discretisation.MeshVariable("cell_size_probe", mesh, 1, degree=1) + return mesh + + +def _cell_size_diagnostics(): + """The three accessors, plus a global sum that no single one of them sees. + + The sum is the sharp one: a per-cell change that leaves the extremes alone + still moves it, and the old field's per-cell values differed while its + minimum happened not to. + """ + mesh = _box() + radii = np.asarray(mesh._cell_radii).reshape(-1) + local_sum = float(radii.sum()) + global_sum = uw.mpi.comm.allreduce(local_sum, op=MPI.SUM) \ + if uw.mpi.size > 1 else local_sum + + values = ( + mesh.get_min_radius(), + mesh.get_max_radius(), + mesh.get_mean_radius(), + global_sum, + ) + return values, mesh_fingerprint(mesh) + + +@pytest.mark.mpi(min_size=2) +def test_cell_size_is_the_same_at_every_rank_count(): + """The accessors and the per-cell sum reproduce their own np=1 answer.""" + values, fingerprint = _cell_size_diagnostics() + compare(values, serial_reference(__file__, "cell_size"), + rtols=RTOLS, labels=LABELS, fingerprint=fingerprint, + what="cell size / radius accessors") + + +@pytest.mark.mpi(min_size=2) +def test_every_rank_agrees_on_the_accessors(): + """The weaker property, kept because it is the one the allreduce gives. + + An allreduce makes an answer identical on every rank; it does not make it + identical at every rank count if the values being reduced are themselves + partition-dependent. Failing this while passing the test above would mean + the reduction is broken rather than its input, so keeping both separates + the two. + """ + mesh = _box() + for name in ("min", "max", "mean"): + got = getattr(mesh, f"get_{name}_radius")() + everyone = uw.mpi.comm.allgather(got) + assert len(set(everyone)) == 1, ( + f"get_{name}_radius() differs between ranks: {everyone}" + ) + + +if __name__ == "__main__": + import sys + + _kind = sys.argv[1] if len(sys.argv) > 1 else "cell_size" + if _kind == "cell_size": + _values, _fingerprint = _cell_size_diagnostics() + emit(_values, _fingerprint) + else: + raise SystemExit(f"unknown kind {_kind!r}") diff --git a/tests/test_0010_cell_size_geometry.py b/tests/test_0010_cell_size_geometry.py index 9e955328e..f7c80bd5c 100644 --- a/tests/test_0010_cell_size_geometry.py +++ b/tests/test_0010_cell_size_geometry.py @@ -1,7 +1,14 @@ -"""Issue #687: cell_size is an own-cell geometric quantity, including after deform. +"""``mesh.cell_size()`` is each cell's ``volume**(1/dim)``, and it tracks deformation. -The independent oracle reads vertex coordinates through the coordinate section; -it does not use the mesh's cached radii or centroid kd-tree. Run serial and MPI. +The oracle here is computed from the vertex coordinates, not from PETSc, so this +is a check and not a restatement of the implementation. For simplices that is +the determinant volume of the cell; for the structured boxes it is the analytic +cell volume, scaled by the determinant of the affine map when the mesh is +deformed. + +Partition independence is NOT tested here -- it cannot be, at one rank count. +``tests/parallel/test_1077`` compares the field cell by cell against its own +serial answer, and ``test_1078`` does the same for the three radius accessors. """ import numpy as np @@ -9,58 +16,87 @@ import underworld3 as uw -pytestmark = [pytest.mark.level_1, pytest.mark.tier_b] +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + +# The deformation applied below, and its Jacobian determinant. Volumes scale by +# |det|, so lengths scale by |det|**(1/dim). +AFFINE = np.array([[1.7, 0.2], [0.0, 1.0]]) +AFFINE_DET = 1.7 + + +def _simplex_volume_from_vertices(mesh): + """Each simplex's volume from its own vertex coordinates. + Triangle: ``|det[v1-v0, v2-v0]| / 2``. Tetrahedron: ``|det[...]| / 6``. + """ + dim = mesh.dim + cell_start, cell_end = mesh.dm.getHeightStratum(0) + point_start, _point_end = mesh.dm.getDepthStratum(0) + + volumes = np.empty(cell_end - cell_start) + factorial = 2.0 if dim == 2 else 6.0 + for cell in range(cell_end - cell_start): + corners = mesh.dm.getTransitiveClosure(cell)[0][-(dim + 1):] + coords = mesh._coords[corners - point_start] + edges = coords[1:] - coords[0] + volumes[cell] = abs(np.linalg.det(edges)) / factorial + return volumes + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_simplex_cell_size_is_the_cube_root_of_its_own_volume(dim): + """Against a determinant volume computed from the cell's vertices.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0,) * dim, maxCoords=(1.0,) * dim, qdegree=3, + cellSize=0.25, regular=False) + mesh.cell_size() + actual = np.asarray(mesh._cell_size_variable.array[:, 0, 0]) + expected = _simplex_volume_from_vertices(mesh) ** (1.0 / dim) -def _vertex_rms(mesh): - dm = mesh.dm - section = dm.getCoordinateDM().getLocalSection() - coordinates = dm.getCoordinatesLocal().array - start, end = dm.getHeightStratum(0) - first_vertex, last_vertex = dm.getDepthStratum(0) - radii = [] - for cell in range(start, end): - vertices = [int(point) for point in dm.getTransitiveClosure(cell)[0] - if first_vertex <= point < last_vertex] - points = np.array([coordinates[section.getOffset(v):section.getOffset(v) + mesh.cdim] - for v in vertices]) - radii.append(np.sqrt(np.mean(np.sum((points - points.mean(axis=0)) ** 2, axis=1)))) - return np.asarray(radii) + error = float(np.abs(actual - expected).max(initial=0.0)) + uw.pprint(f"CELL_SIZE_GEOMETRY simplex dim={dim} max_error={error:.12g}") + assert error < 1.0e-12, error @pytest.mark.parametrize("dim", [2, 3]) -@pytest.mark.parametrize("simplex", [True, False], ids=["simplex", "tensor"]) -def test_cell_size_matches_own_vertices_and_tracks_deform(dim, simplex): - geometry = dict(minCoords=(0.0,) * dim, maxCoords=(1.0,) * dim, qdegree=3) - mesh = (uw.meshing.UnstructuredSimplexBox(**geometry, cellSize=0.25, regular=False) - if simplex else uw.meshing.StructuredQuadBox(**geometry, elementRes=(4,) * dim)) +def test_structured_cell_size_is_the_analytic_value_and_tracks_deformation(dim): + """A regular box has an exact answer, and an affine deform scales it. + + ``elementRes=4`` on the unit box gives cells of side 0.25, so + ``volume**(1/dim)`` is 0.25 whatever the dimension. Applying a linear map + multiplies every cell volume by ``|det|``, hence every length by + ``|det|**(1/dim)`` -- which the field must follow after ``deform``. + """ + mesh = uw.meshing.StructuredQuadBox( + minCoords=(0.0,) * dim, maxCoords=(1.0,) * dim, qdegree=3, + elementRes=(4,) * dim) mesh.cell_size() field = mesh._cell_size_variable - errors = [] - for phase in ("initial", "deformed"): - if phase == "deformed": - coordinates = np.array(mesh.X.coords) - coordinates[:, 0] = 1.7 * coordinates[:, 0] + 0.2 * coordinates[:, 1] - mesh.deform(coordinates) - expected = _vertex_rms(mesh) - actual = np.asarray(field.array[:, 0, 0]) - shapes_match = actual.shape == expected.shape - assert all(uw.mpi.comm.allgather(shapes_match)), (actual.shape, expected.shape) - local_error = float(np.abs(actual - expected).max(initial=0.0)) - error = max(uw.mpi.comm.allgather(local_error)) - errors.append(error) - uw.pprint(f"CELL_SIZE_GEOMETRY dim={dim} simplex={simplex} phase={phase} " - f"ranks={uw.mpi.size} max_error={error:.12g}") - assert max(errors) < 1e-12, errors - - -def test_regular_square_cell_size_keeps_global_radius(): + + actual = np.asarray(field.array[:, 0, 0]) + error = float(np.abs(actual - 0.25).max(initial=0.0)) + assert max(uw.mpi.comm.allgather(error)) < 1.0e-12, error + + coordinates = np.array(mesh.X.coords) + coordinates[:, 0] = AFFINE[0, 0] * coordinates[:, 0] + AFFINE[0, 1] * coordinates[:, 1] + mesh.deform(coordinates) + + expected = 0.25 * AFFINE_DET ** (1.0 / dim) + actual = np.asarray(field.array[:, 0, 0]) + error = float(np.abs(actual - expected).max(initial=0.0)) + uw.pprint(f"CELL_SIZE_GEOMETRY tensor dim={dim} deformed max_error={error:.12g}") + assert max(uw.mpi.comm.allgather(error)) < 1.0e-12, (error, expected) + + +def test_the_global_minimum_agrees_with_the_field(): + """``get_min_radius()`` reduces the same field ``cell_size()`` exposes. + + On a regular box every cell is the same size, so the global minimum is that + size -- which also pins the reduction to the analytic value rather than to + whatever the field happens to hold. + """ mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4), qdegree=2) - legacy = np.array(mesh._radii) - global_radius = mesh.get_min_radius() mesh.cell_size() - expected = np.sqrt(2.0) / 8.0 - error = float(np.abs(np.asarray(mesh._cell_size_variable.array) - expected).max(initial=0.0)) - assert max(uw.mpi.comm.allgather(error)) < 1e-12 - assert global_radius == pytest.approx(expected, rel=1e-12) - assert all(uw.mpi.comm.allgather(np.array_equal(mesh._radii, legacy))) + assert mesh.get_min_radius() == pytest.approx(0.25, rel=1.0e-12) + assert float(np.asarray(mesh._cell_size_variable.array).max()) == pytest.approx( + 0.25, rel=1.0e-12) diff --git a/tests/test_1060_nitsche_freeslip.py b/tests/test_1060_nitsche_freeslip.py index fe436d00d..4089f2875 100644 --- a/tests/test_1060_nitsche_freeslip.py +++ b/tests/test_1060_nitsche_freeslip.py @@ -67,8 +67,8 @@ def _solve_freeslip_box(method, res=8): stokes.add_natural_bc(1e4 * Gamma.dot(v.sym) * Gamma, "Top") stokes.add_natural_bc(1e4 * Gamma.dot(v.sym) * Gamma, "Bottom") elif method == "nitsche": - stokes.add_nitsche_bc(0.0, "Top", gamma=10.0) - stokes.add_nitsche_bc(0.0, "Bottom", gamma=10.0) + stokes.add_nitsche_bc(0.0, "Top") + stokes.add_nitsche_bc(0.0, "Bottom") else: raise ValueError(f"Unknown method: {method}") @@ -147,4 +147,4 @@ def test_nitsche_better_than_penalty_constraint(self, solutions): max_vn_pen = np.max(np.abs(v_pen[top_pen, 1])) if np.any(top_pen) else 0 print(f"Normal velocity on top: Nitsche={max_vn_nit:.4e}, Penalty={max_vn_pen:.4e}") - # Nitsche at gamma=10 should be comparable or better than penalty at 1e4 + # Nitsche at the default gamma should be comparable or better than penalty at 1e4 diff --git a/tests/test_1065_nitsche_local_h.py b/tests/test_1065_nitsche_local_h.py index e6c225611..b643b46b4 100644 --- a/tests/test_1065_nitsche_local_h.py +++ b/tests/test_1065_nitsche_local_h.py @@ -152,7 +152,7 @@ def _box_wobble(X0, amp): # -------------------------------------------------------------------------- def test_cell_size_is_local_per_cell(): """``mesh.cell_size()`` is a per-cell field equal to each cell's - own-vertex RMS size (``mesh._cell_radii``), not the single global minimum.""" + per-cell size (``mesh._cell_radii``, PETSc's ``volume**(1/dim)``), not the single global minimum.""" mesh = _graded_box() h = mesh.cell_size() # sympy symbol -> backed by a P0 field field = np.asarray(mesh._cell_size_variable.array[:, 0, 0]).reshape(-1) @@ -225,7 +225,7 @@ def test_cell_size_tracks_deformation(): # -------------------------------------------------------------------------- # 3. local-h still solves free-slip correctly (back-compat / correctness) # -------------------------------------------------------------------------- -def _solve_freeslip(mesh, method, gamma=10.0): +def _solve_freeslip(mesh, method, gamma=12.5): v = uw.discretisation.MeshVariable( "U", mesh, mesh.dim, degree=2, vtype=uw.VarType.VECTOR) p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1)