Output metadata in serial so VLEN strings can be used

gprMax/output_controllers/geometry_view_lines.py

@@ -181,5 +181,18 @@ class MPIGeometryViewLines(GeometryViewLines[MPIGrid]):
             self.cell_offsets,
             comm=self.grid_view.comm,
         ) as f:
-            self.metadata.write_to_vtkhdf(f)
             f.add_cell_data("Material", self.material_data, self.grid_view.offset)
+
+        # Write metadata in serial as it contains variable length
+        # strings which currently cannot be written by HDF5 using
+        # parallel I/O
+        if self.grid_view.comm.rank == 0:
+            with VtkUnstructuredGrid(
+                self.filename,
+                self.points,
+                self.cell_types,
+                self.connectivity,
+                self.cell_offsets,
+                mode="r+",
+            ) as f:
+                self.metadata.write_to_vtkhdf(f)

gprMax/output_controllers/geometry_view_voxels.py

@@ -98,5 +98,13 @@ class MPIGeometryViewVoxels(GeometryViewVoxels[MPIGrid]):
             self.spacing,
             comm=self.grid_view.comm,
         ) as f:
-            self.metadata.write_to_vtkhdf(f)
             f.add_cell_data("Material", self.material_data, self.grid_view.offset)
+
+        # Write metadata in serial as it contains variable length
+        # strings which currently cannot be written by HDF5 using
+        # parallel I/O
+        if self.grid_view.comm.rank == 0:
+            with VtkImageData(
+                self.filename, self.grid_view.global_size, self.origin, self.spacing, mode="r+"
+            ) as f:
+                self.metadata.write_to_vtkhdf(f)

gprMax/output_controllers/geometry_views.py

@@ -178,7 +178,7 @@ class Metadata(Generic[GridType]):
         file_handler.add_field_data("dx_dy_dz", self.dx_dy_dz)
         file_handler.add_field_data("nx_ny_nz", self.nx_ny_nz)
 
-        self.write_material_ids(file_handler)
+        file_handler.add_field_data("material_ids", self.materials)
 
         if not self.materials_only:
             if self.pml_thickness is not None:
@@ -192,9 +192,6 @@ class Metadata(Generic[GridType]):
                 file_handler.add_field_data("receiver_ids", self.receiver_ids)
                 file_handler.add_field_data("receivers", self.receiver_positions)
 
-    def write_material_ids(self, file_handler: VtkHdfFile):
-        file_handler.add_field_data("material_ids", self.materials)
-
     def pml_gv_comment(self) -> Optional[npt.NDArray[np.int64]]:
         grid = self.grid
 
@@ -291,22 +288,3 @@ class MPIMetadata(Metadata[MPIGrid]):
         objects = dict(sorted(objects.items()))
 
         return (list(objects.keys()), np.array(list(objects.values()))) if objects else None
-
-    def write_material_ids(self, file_handler: VtkHdfFile):
-        assert isinstance(self.grid_view, MPIGridView)
-
-        # Only rank 0 has all the material data. However, creating the
-        # 'material_ids' dataset is a collective operation, so all ranks
-        # need to know the shape and datatype of the dataset.
-        if self.materials is None:
-            buffer = np.empty(2, dtype=np.int32)
-        else:
-            shape = len(self.materials)
-            max_length = max([len(m) for m in self.materials])
-            buffer = np.array([shape, max_length], dtype=np.int32)
-
-        self.grid_view.comm.Bcast([buffer, MPI.INT32_T])
-        shape, max_length = buffer
-        dtype = h5py.string_dtype(length=int(max_length))
-
-        file_handler.add_field_data("material_ids", self.materials, shape=(shape,), dtype=dtype)