Swarm/mesh checkpoint quiescence, .msh coordinate-system default, rank-symmetric stats temporaries

src/underworld3/discretisation/discretisation_mesh.py

@@ -795,7 +795,22 @@ def _load_dm_from_file(
                     coord_type_dict = json.loads(json_str)
                     coordinate_system_type = CoordinateSystemType(coord_type_dict["value"])
                 except KeyError:
-                    pass
+                    # A checkpointed mesh normally carries its coordinate
+                    # system; a missing entry means the metadata was lost
+                    # (or the file predates it). The construction defaults
+                    # to CARTESIAN (#397) — say so, because for an annulus
+                    # or spherical checkpoint that label is wrong and would
+                    # otherwise propagate silently into re-written files.
+                    if coordinate_system_type is None and uw.mpi.rank == 0:
+                        import warnings
+
+                        warnings.warn(
+                            f"{plex_or_meshfile} has no coordinate_system_type "
+                            "metadata; defaulting to CARTESIAN. Pass "
+                            "coordinate_system_type= explicitly if this mesh "
+                            "is not Cartesian.",
+                            stacklevel=2,
+                        )
 
                 regions = None
                 try:
@@ -1223,6 +1238,14 @@ def _setup_symbolic_coordinates(self, coordinate_system_type):
         # A unique set of vectors / names for each mesh instance
         #
 
+        # A mesh constructed without an explicit coordinate system (e.g.
+        # loaded directly from a .msh file, or from an h5 checkpoint with no
+        # coordinate metadata) is Cartesian — the default every uw.meshing
+        # constructor passes. Leaving None here crashed mesh.write()'s
+        # metadata block (issue #397).
+        if coordinate_system_type is None:
+            coordinate_system_type = CoordinateSystemType.CARTESIAN
+
         self.CoordinateSystemType = coordinate_system_type
 
         from sympy.vector import CoordSys3D
@@ -4815,58 +4838,63 @@ def write(
         # Use preferred selective_ranks pattern for metadata operations
         with uw.selective_ranks(0) as should_execute:
             if should_execute:
-                f = h5py.File(filename, "a")
-                g = f.create_group("metadata")
-
-                boundaries_dict = {i.name: i.value for i in self.boundaries}
-                g.attrs["boundaries"] = json.dumps(boundaries_dict)
-
-                if self.regions is not None:
-                    regions_dict = {i.name: i.value for i in self.regions}
-                    g.attrs["regions"] = json.dumps(regions_dict)
-
-                coordinates_type_dict = {
-                    "name": self.CoordinateSystemType.name,
-                    "value": self.CoordinateSystemType.value,
-                }
-                g.attrs["coordinate_system_type"] = json.dumps(coordinates_type_dict)
-
-                # Save ellipsoid metadata for geographic meshes
-                if hasattr(self.CoordinateSystem, "ellipsoid"):
-                    ellipsoid_ser = {}
-                    for k, v in self.CoordinateSystem.ellipsoid.items():
-                        if hasattr(v, "to"):  # uw.quantity
-                            ellipsoid_ser[k] = {
-                                "value": float(v.magnitude),
-                                "unit": str(v.units),
-                            }
-                        else:
-                            ellipsoid_ser[k] = v
-                    g.attrs["ellipsoid"] = json.dumps(ellipsoid_ser)
-
-                # Add coordinate units metadata
-                if hasattr(self, "coordinate_units"):
-                    coord_units_dict = {
-                        "coordinate_units": str(self.coordinate_units),
-                        "coordinate_dimensionality": (
-                            str(self.coordinate_dimensionality)
-                            if hasattr(self, "coordinate_dimensionality")
-                            else None
-                        ),
-                        "length_scale": (
-                            str(self.length_scale) if hasattr(self, "length_scale") else None
-                        ),
-                        "mesh_type": type(self).__name__,
-                        "dimension": self.dim,
-                    }
-                    g.attrs["coordinate_units"] = json.dumps(coord_units_dict)
+                # Context manager: an exception mid-block must not leak the
+                # handle (a live handle keeps the HDF5 lock held).
+                with h5py.File(filename, "a") as f:
+                    g = f.create_group("metadata")
 
-                # Number of coarse multigrid levels in the hierarchy (= number
-                # of refinements from the stored coarsest level up to the fine
-                # mesh). Used on reload to rebuild the intermediate levels.
-                g.attrs["hierarchy_coarse_levels"] = len(self.dm_hierarchy) - 1
+                    boundaries_dict = {i.name: i.value for i in self.boundaries}
+                    g.attrs["boundaries"] = json.dumps(boundaries_dict)
 
-                f.close()
+                    if self.regions is not None:
+                        regions_dict = {i.name: i.value for i in self.regions}
+                        g.attrs["regions"] = json.dumps(regions_dict)
+
+                    coordinates_type_dict = {
+                        "name": self.CoordinateSystemType.name,
+                        "value": self.CoordinateSystemType.value,
+                    }
+                    g.attrs["coordinate_system_type"] = json.dumps(coordinates_type_dict)
+
+                    # Save ellipsoid metadata for geographic meshes
+                    if hasattr(self.CoordinateSystem, "ellipsoid"):
+                        ellipsoid_ser = {}
+                        for k, v in self.CoordinateSystem.ellipsoid.items():
+                            if hasattr(v, "to"):  # uw.quantity
+                                ellipsoid_ser[k] = {
+                                    "value": float(v.magnitude),
+                                    "unit": str(v.units),
+                                }
+                            else:
+                                ellipsoid_ser[k] = v
+                        g.attrs["ellipsoid"] = json.dumps(ellipsoid_ser)
+
+                    # Add coordinate units metadata
+                    if hasattr(self, "coordinate_units"):
+                        coord_units_dict = {
+                            "coordinate_units": str(self.coordinate_units),
+                            "coordinate_dimensionality": (
+                                str(self.coordinate_dimensionality)
+                                if hasattr(self, "coordinate_dimensionality")
+                                else None
+                            ),
+                            "length_scale": (
+                                str(self.length_scale) if hasattr(self, "length_scale") else None
+                            ),
+                            "mesh_type": type(self).__name__,
+                            "dimension": self.dim,
+                        }
+                        g.attrs["coordinate_units"] = json.dumps(coord_units_dict)
+
+                    # Number of coarse multigrid levels in the hierarchy (= number
+                    # of refinements from the stored coarsest level up to the fine
+                    # mesh). Used on reload to rebuild the intermediate levels.
+                    g.attrs["hierarchy_coarse_levels"] = len(self.dm_hierarchy) - 1
+
+        # Same quiescence contract as Swarm.save (issue #330): every rank
+        # waits for rank 0's metadata append, so an immediate reopen of the
+        # mesh file cannot hit HDF5 file locking.
+        uw.mpi.barrier()
 
         # Persist the geometric-multigrid (FMG) hierarchy as a SINGLE sidecar
         # holding the coarsest level only. On reload the intermediate coarse

src/underworld3/discretisation/discretisation_mesh_variables.py

@@ -986,28 +986,33 @@ def save(
         # Use preferred selective_ranks pattern for unit metadata
         with uw.selective_ranks(0) as should_execute:
             if should_execute:
-                f = h5py.File(filename, "a")
-
-                # Create or get metadata group
-                if "metadata" not in f:
-                    g = f.create_group("metadata")
-                else:
-                    g = f["metadata"]
-
-                # Add variable unit metadata
-                var_metadata = {
-                    "units": str(self.units) if hasattr(self, "units") and self.units else None,
-                    "dimensionality": (
-                        str(self.dimensionality) if hasattr(self, "dimensionality") else None
-                    ),
-                    "units_backend": "pint" if self.has_units else None,
-                    "num_components": self.num_components,
-                    "variable_type": str(self.vtype),
-                    "variable_name": self.name,
-                }
-
-                g.attrs[f"variable_{self.clean_name}_units"] = json.dumps(var_metadata)
-                f.close()
+                # Context manager: an exception mid-block must not leak the
+                # handle (a live handle keeps the HDF5 lock held).
+                with h5py.File(filename, "a") as f:
+                    # Create or get metadata group
+                    if "metadata" not in f:
+                        g = f.create_group("metadata")
+                    else:
+                        g = f["metadata"]
+
+                    # Add variable unit metadata
+                    var_metadata = {
+                        "units": str(self.units) if hasattr(self, "units") and self.units else None,
+                        "dimensionality": (
+                            str(self.dimensionality) if hasattr(self, "dimensionality") else None
+                        ),
+                        "units_backend": "pint" if self.has_units else None,
+                        "num_components": self.num_components,
+                        "variable_type": str(self.vtype),
+                        "variable_name": self.name,
+                    }
+
+                    g.attrs[f"variable_{self.clean_name}_units"] = json.dumps(var_metadata)
+
+        # Same quiescence contract as Swarm.save (issue #330): every rank
+        # waits for rank 0's metadata append, so an immediate reopen
+        # (read_timestep after write_timestep) cannot hit HDF5 file locking.
+        uw.mpi.barrier()
 
         lvec = self.mesh.dm.getCoordinates()
 
@@ -3126,8 +3131,20 @@ def _vector_stats(self):
         """Statistics for vector variables using magnitude."""
         import numpy as np
 
-        # Create temporary scalar variable for magnitude
-        magnitude_var = _BaseMeshVariable(f"_temp_mag_{id(self)}", self.mesh, 1, degree=self.degree)
+        # Temporary scalar variable for the magnitude. The name suffix must
+        # be rank-symmetric: variable creation/destruction performs
+        # collective DM operations keyed by field name, and id(self) is a
+        # rank-local address that differs across ranks (issue #384).
+        temp_name = f"_temp_mag_{self.clean_name}"
+        # A name collision would silently ALIAS an existing variable
+        # (creation returns the registered object), and the stats pass
+        # would then overwrite its data and deregister it — refuse instead.
+        if temp_name in self.mesh.vars:
+            raise RuntimeError(
+                f"Cannot compute vector stats: a variable named '{temp_name}' "
+                "already exists on this mesh (reserved as a stats temporary)."
+            )
+        magnitude_var = _BaseMeshVariable(temp_name, self.mesh, 1, degree=self.degree)
 
         try:
             # Compute magnitude: |v| = sqrt(v·v)
@@ -3164,9 +3181,17 @@ def _tensor_stats(self):
         """Statistics for tensor variables using Frobenius norm."""
         import numpy as np
 
-        # Create temporary scalar variable for Frobenius norm
+        # Temporary scalar variable for the Frobenius norm — rank-symmetric
+        # name suffix and collision refusal for the same reasons as
+        # _vector_stats (issue #384).
+        temp_name = f"_temp_frob_{self.clean_name}"
+        if temp_name in self.mesh.vars:
+            raise RuntimeError(
+                f"Cannot compute tensor stats: a variable named '{temp_name}' "
+                "already exists on this mesh (reserved as a stats temporary)."
+            )
         frobenius_var = uw.discretisation.MeshVariable(
-            f"_temp_frob_{id(self)}", self.mesh, 1, degree=self.degree
+            temp_name, self.mesh, 1, degree=self.degree
         )
 
         try:

src/underworld3/swarm.py

@@ -2016,6 +2016,11 @@ def save(
                     if "data" in h5f:
                         h5f["data"].attrs["units_metadata"] = json.dumps(swarm_metadata)
 
+        # Same quiescence contract as Swarm.save (issue #330): all ranks
+        # wait for rank 0's metadata append so an immediate reopen cannot
+        # hit HDF5 file locking.
+        comm.barrier()
+
         return
 
     @timing.routine_timer_decorator
@@ -4087,11 +4092,6 @@ def save(
             del points_data_copy
 
         ## Add swarm coordinate unit metadata to the file
-        # TODO(BUG): save() returns on non-zero ranks while rank 0 is still
-        # appending this metadata; a rank that immediately reopens the file
-        # (e.g. read_timestep straight after write_timestep) hits HDF5 file
-        # locking (BlockingIOError, errno 35). A trailing barrier would make
-        # the file quiescent on return. Found while reproducing issue #324.
         import json
 
         # Use preferred selective_ranks pattern for coordinate metadata
@@ -4121,6 +4121,13 @@ def save(
                             swarm_coord_metadata
                         )
 
+        # The file must be quiescent when save() returns on EVERY rank:
+        # without this barrier, non-zero ranks return while rank 0 still
+        # holds the file open for the metadata append, and an immediate
+        # reopen (e.g. read_timestep right after write_timestep) hits HDF5
+        # file locking (BlockingIOError, errno 35) — issue #330.
+        comm.barrier()
+
         return
 
     @timing.routine_timer_decorator

tests/test_0003_save_load.py

@@ -277,7 +277,10 @@ def test_write_timestep_missing_directory_raises(tmp_path):
     """
     import underworld3 as uw
 
-    mesh = uw.meshing.StructuredQuadBox(elementRes=(2, 2))
+    # 4x4 keeps every rank populated up to np4 — a 2x2 mesh starves ranks
+    # of cells and crashes in mesh construction (issue #399), which then
+    # desynchronises every test after this one in an MPI run.
+    mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4))
     missing_dir = tmp_path / "does_not_exist"
 
     with pytest.raises(RuntimeError, match="does not exist"):
@@ -372,3 +375,24 @@ def test_read_timestep_refuses_frame_mismatched_file(tmp_path):
     w = uw.discretisation.MeshVariable("phi269b", mesh, 1, degree=0, continuous=False)
     with pytest.raises(RuntimeError, match="coordinate frames do not match"):
         w.read_timestep("chk269g", "phi269", 0, outputPath=tmp_path)
+
+
+def test_mesh_from_msh_file_writes(tmp_path):
+    """#397: a mesh loaded directly from a .msh file must be writable.
+
+    The coordinate system defaults to CARTESIAN (matching every uw.meshing
+    constructor) instead of leaving ``CoordinateSystemType=None`` and
+    crashing ``mesh.write()``'s metadata block.
+    """
+    import underworld3 as uw
+
+    tmp_path = _shared_tmp_path(tmp_path, uw)
+
+    msh = f"{tmp_path}/box397.msh"
+    scratch = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25, filename=msh)
+    del scratch
+
+    mesh = uw.discretisation.Mesh(msh)
+    assert mesh.CoordinateSystemType is uw.coordinates.CoordinateSystemType.CARTESIAN
+    mesh.write(f"{tmp_path}/box397.h5")