diff --git a/src/underworld3/discretisation/discretisation_mesh.py b/src/underworld3/discretisation/discretisation_mesh.py index 20c8516b..7cc36e9b 100644 --- a/src/underworld3/discretisation/discretisation_mesh.py +++ b/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 diff --git a/src/underworld3/discretisation/discretisation_mesh_variables.py b/src/underworld3/discretisation/discretisation_mesh_variables.py index 02f4ae8b..188fca2a 100644 --- a/src/underworld3/discretisation/discretisation_mesh_variables.py +++ b/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: diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py index 90f912bc..06c52276 100644 --- a/src/underworld3/swarm.py +++ b/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 diff --git a/tests/test_0003_save_load.py b/tests/test_0003_save_load.py index 1b91b743..0b2d1839 100644 --- a/tests/test_0003_save_load.py +++ b/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")