Update geometry views for parallel build

- Sources and receivers now set their ID using a static point so the
  value is the same when using the MPI implementation.
- The model has new functions to create geometry objects and geometry
  views.
- Grid views create their own MPI communicator. This means when creating
  a grid view, or an object containing a grid view, it is the creators
  responsibility to call `comm.Split(MPI.UNDEFINED)` for ranks not
  creating the grid view.
- There is a new MPIMetadata class for managing geometry view metadata.
- Ranks can create field data datasets in VTKHDF files without writing
  any data.
- The VtkHdfFile class now correctly handles string datatypes.

gprMax/model.py

@@ -30,6 +30,9 @@ from colorama import Fore, Style, init
 from gprMax.grid.cuda_grid import CUDAGrid
 from gprMax.grid.opencl_grid import OpenCLGrid
 from gprMax.output_controllers.geometry_objects import GeometryObject
+from gprMax.output_controllers.geometry_view_lines import GeometryViewLines
+from gprMax.output_controllers.geometry_view_voxels import GeometryViewVoxels
+from gprMax.output_controllers.geometry_views import GeometryView, save_geometry_views
 from gprMax.subgrids.grid import SubGridBaseGrid
 
 init()
@@ -39,7 +42,6 @@ from tqdm import tqdm
 import gprMax.config as config
 
 from .fields_outputs import write_hdf5_outputfile
-from .geometry_outputs import GeometryView, save_geometry_views
 from .grid.fdtd_grid import FDTDGrid
 from .snapshots import save_snapshots
 from .utilities.host_info import mem_check_build_all, mem_check_run_all, set_omp_threads
@@ -172,6 +174,11 @@ class Model:
         return grid
 
     def set_size(self, size: npt.NDArray[np.int32]):
+        """Set size of the model.
+
+        Args:
+            size: Array to set the size (3 dimensional).
+        """
         self.nx, self.ny, self.nz = size
 
     def add_geometry_object(
@@ -181,12 +188,90 @@ class Model:
         stop: npt.NDArray[np.int32],
         basefilename: str,
     ) -> Optional[GeometryObject]:
+        """Add a geometry object to the model.
+
+        Args:
+            grid: Grid to create a geometry object for.
+            start: Lower extent of the geometry object (3 dimensional).
+            stop: Upper extent of the geometry object (3 dimensional).
+            basefilename: Output filename of the geometry object.
+
+        Returns:
+            geometry_object: The created geometry object.
+        """
         geometry_object = GeometryObject(
             grid, start[0], start[1], start[2], stop[0], stop[1], stop[2], basefilename
         )
         self.geometryobjects.append(geometry_object)
         return geometry_object
 
+    def add_geometry_view_voxels(
+        self,
+        grid: FDTDGrid,
+        start: npt.NDArray[np.int32],
+        stop: npt.NDArray[np.int32],
+        dl: npt.NDArray[np.int32],
+        filename: str,
+    ) -> Optional[GeometryViewVoxels]:
+        """Add a voxel geometry view to the model.
+
+        Args:
+            grid: Grid to create a geometry view for.
+            start: Lower extent of the geometry view (3 dimensional).
+            stop: Upper extent of the geometry view (3 dimensional).
+            dl: Discritisation of the geometry view (3 dimensional).
+            filename: Output filename of the geometry view.
+
+        Returns:
+            geometry_view: The created geometry view.
+        """
+        geometry_view = GeometryViewVoxels(
+            start[0],
+            start[1],
+            start[2],
+            stop[0],
+            stop[1],
+            stop[2],
+            dl[0],
+            dl[1],
+            dl[2],
+            filename,
+            grid,
+        )
+        self.geometryviews.append(geometry_view)
+        return geometry_view
+
+    def add_geometry_view_lines(
+        self,
+        grid: FDTDGrid,
+        start: npt.NDArray[np.int32],
+        stop: npt.NDArray[np.int32],
+        filename: str,
+    ) -> Optional[GeometryViewLines]:
+        """Add a lines geometry view to the model.
+
+        Args:
+            grid: Grid to create a geometry view for.
+            start: Lower extent of the geometry view (3 dimensional).
+            stop: Upper extent of the geometry view (3 dimensional).
+            filename: Output filename of the geometry view.
+
+        Returns:
+            geometry_view: The created geometry view.
+        """
+        geometry_view = GeometryViewLines(
+            start[0],
+            start[1],
+            start[2],
+            stop[0],
+            stop[1],
+            stop[2],
+            filename,
+            grid,
+        )
+        self.geometryviews.append(geometry_view)
+        return geometry_view
+
     def build(self):
         """Builds the Yee cells for a model."""
 

gprMax/mpi_model.py

@@ -10,6 +10,8 @@ from gprMax.fields_outputs import write_hdf5_outputfile
 from gprMax.grid.mpi_grid import MPIGrid
 from gprMax.model import Model
 from gprMax.output_controllers.geometry_objects import MPIGeometryObject
+from gprMax.output_controllers.geometry_view_lines import MPIGeometryViewLines
+from gprMax.output_controllers.geometry_view_voxels import MPIGeometryViewVoxels
 from gprMax.snapshots import save_snapshots
 
 logger = logging.getLogger(__name__)
@@ -67,16 +69,116 @@ class MPIModel(Model):
         stop: npt.NDArray[np.int32],
         basefilename: str,
     ) -> Optional[MPIGeometryObject]:
-        comm = grid.create_sub_communicator(start, stop)
+        """Add a geometry object to the model.
 
-        if comm is None:
-            return None
-        else:
+        Args:
+            grid: Grid to create a geometry object for
+            start: Lower extent of the geometry object (3 dimensional)
+            stop: Upper extent of the geometry object (3 dimensional)
+            basefilename: Output filename of the geometry object
+
+        Returns:
+            geometry_object: The new geometry object or None if no
+                geometry object was created.
+        """
+        if grid.local_bounds_overlap_grid(start, stop):
             geometry_object = MPIGeometryObject(
-                grid, start[0], start[1], start[2], stop[0], stop[1], stop[2], basefilename, comm
+                grid, start[0], start[1], start[2], stop[0], stop[1], stop[2], basefilename
             )
             self.geometryobjects.append(geometry_object)
             return geometry_object
+        else:
+            # The MPIGridView created by the MPIGeometryObject will
+            # create a new communicator using MPI_Split. Calling this
+            # here prevents deadlock if not all ranks create the new
+            # MPIGeometryObject.
+            grid.comm.Split(MPI.UNDEFINED)
+            return None
+
+    def add_geometry_view_voxels(
+        self,
+        grid: MPIGrid,
+        start: npt.NDArray[np.int32],
+        stop: npt.NDArray[np.int32],
+        dl: npt.NDArray[np.int32],
+        filename: str,
+    ) -> Optional[MPIGeometryViewVoxels]:
+        """Add a voxel geometry view to the model.
+
+        Args:
+            grid: Grid to create a geometry view for.
+            start: Lower extent of the geometry view (3 dimensional).
+            stop: Upper extent of the geometry view (3 dimensional).
+            dl: Discritisation of the geometry view (3 dimensional).
+            filename: Output filename of the geometry view.
+
+        Returns:
+            geometry_view: The new geometry view or None if no geometry
+                view was created.
+        """
+        if grid.local_bounds_overlap_grid(start, stop):
+            geometry_view = MPIGeometryViewVoxels(
+                start[0],
+                start[1],
+                start[2],
+                stop[0],
+                stop[1],
+                stop[2],
+                dl[0],
+                dl[1],
+                dl[2],
+                filename,
+                grid,
+            )
+            self.geometryviews.append(geometry_view)
+            return geometry_view
+        else:
+            # The MPIGridView created by MPIGeometryViewVoxels will
+            # create a new communicator using MPI_Split. Calling this
+            # here prevents deadlock if not all ranks create the new
+            # MPIGeometryViewVoxels.
+            grid.comm.Split(MPI.UNDEFINED)
+            return None
+
+    def add_geometry_view_lines(
+        self,
+        grid: MPIGrid,
+        start: npt.NDArray[np.int32],
+        stop: npt.NDArray[np.int32],
+        filename: str,
+    ) -> Optional[MPIGeometryViewLines]:
+        """Add a lines geometry view to the model.
+
+        Args:
+            grid: Grid to create a geometry view for.
+            start: Lower extent of the geometry view (3 dimensional).
+            stop: Upper extent of the geometry view (3 dimensional).
+            filename: Output filename of the geometry view.
+
+        Returns:
+            geometry_view: The new geometry view or None if no geometry
+                view was created.
+        """
+        if grid.local_bounds_overlap_grid(start, stop):
+            geometry_view = MPIGeometryViewLines(
+                start[0],
+                start[1],
+                start[2],
+                stop[0],
+                stop[1],
+                stop[2],
+                filename,
+                grid,
+            )
+            self.geometryviews.append(geometry_view)
+            return geometry_view
+        else:
+            # The MPIGridView created by MPIGeometryViewLines will
+            # create a new communicator using MPI_Split. Calling this
+            # here prevents deadlock if not all ranks create the new
+            # MPIGeometryViewLines.
+            grid.comm.Split(MPI.UNDEFINED)
+            return None
 
     def write_output_data(self):
         """Writes output data, i.e. field data for receivers and snapshots to

gprMax/output_controllers/geometry_objects.py

@@ -1,5 +1,6 @@
 from io import TextIOWrapper
 from pathlib import Path
+from typing import Generic
 
 import h5py
 import numpy as np
@@ -11,10 +12,10 @@ from gprMax._version import __version__
 from gprMax.grid.fdtd_grid import FDTDGrid
 from gprMax.grid.mpi_grid import MPIGrid
 from gprMax.materials import Material
-from gprMax.output_controllers.grid_view import GridView, MPIGridView
+from gprMax.output_controllers.grid_view import GridType, GridView, MPIGridView
 
 
-class GeometryObject:
+class GeometryObject(Generic[GridType]):
     """Geometry objects to be written to file."""
 
     @property
@@ -22,7 +23,7 @@ class GeometryObject:
         return GridView
 
     def __init__(
-        self, grid: FDTDGrid, xs: int, ys: int, zs: int, xf: int, yf: int, zf: int, filename: str
+        self, grid: GridType, xs: int, ys: int, zs: int, xf: int, yf: int, zf: int, filename: str
     ):
         """
         Args:
@@ -45,7 +46,7 @@ class GeometryObject:
         self.datawritesize = self.solidsize + self.rigidsize + self.IDsize
 
     @property
-    def grid(self) -> FDTDGrid:
+    def grid(self) -> GridType:
         return self.grid_view.grid
 
     def write_metadata(self, file_handler: h5py.File, title: str):
@@ -115,26 +116,11 @@ class GeometryObject:
                 self.output_material(material, fmaterials)
 
 
-class MPIGeometryObject(GeometryObject):
+class MPIGeometryObject(GeometryObject[MPIGrid]):
     @property
     def GRID_VIEW_TYPE(self) -> type[MPIGridView]:
         return MPIGridView
 
-    def __init__(
-        self,
-        grid: MPIGrid,
-        xs: int,
-        xf: int,
-        ys: int,
-        yf: int,
-        zs: int,
-        zf: int,
-        filename: str,
-        comm: MPI.Cartcomm,
-    ):
-        super().__init__(grid, xs, xf, ys, yf, zs, zf, filename)
-        self.comm = comm
-
     def write_hdf5(self, title: str, pbar: tqdm):
         """Writes a geometry objects file in HDF5 format.
 
@@ -144,7 +130,7 @@ class MPIGeometryObject(GeometryObject):
         """
         assert isinstance(self.grid_view, self.GRID_VIEW_TYPE)
 
-        self.grid_view.initialise_materials(self.comm)
+        self.grid_view.initialise_materials()
 
         ID = self.grid_view.get_ID()
         data = self.grid_view.get_solid().astype(np.int16)
@@ -154,14 +140,10 @@ class MPIGeometryObject(GeometryObject):
         ID = self.grid_view.map_to_view_materials(ID)
         data = self.grid_view.map_to_view_materials(data)
 
-        with h5py.File(self.filename_hdf5, "w", driver="mpio", comm=self.comm) as fdata:
+        with h5py.File(self.filename_hdf5, "w", driver="mpio", comm=self.grid_view.comm) as fdata:
             self.write_metadata(fdata, title)
 
-            dset_slice = (
-                self.grid_view.get_output_slice(0),
-                self.grid_view.get_output_slice(1),
-                self.grid_view.get_output_slice(2),
-            )
+            dset_slice = self.grid_view.get_3d_output_slice()
 
             dset = fdata.create_dataset("/data", self.grid_view.global_size, dtype=data.dtype)
             dset[dset_slice] = data
@@ -178,11 +160,7 @@ class MPIGeometryObject(GeometryObject):
             rigid_H_dataset[:, dset_slice[0], dset_slice[1], dset_slice[2]] = rigidH
             pbar.update(self.rigidsize)
 
-            dset_slice = (
-                self.grid_view.get_output_slice(0, upper_bound_exclusive=False),
-                self.grid_view.get_output_slice(1, upper_bound_exclusive=False),
-                self.grid_view.get_output_slice(2, upper_bound_exclusive=False),
-            )
+            dset_slice = self.grid_view.get_3d_output_slice(upper_bound_exclusive=False)
 
             dset = fdata.create_dataset(
                 "/ID", (6, *(self.grid_view.global_size + 1)), dtype=ID.dtype

gprMax/output_controllers/geometry_view_lines.py

@@ -19,13 +19,13 @@
 import logging
 
 import numpy as np
-from mpi4py import MPI
 
 from gprMax._version import __version__
 from gprMax.cython.geometry_outputs import get_line_properties
+from gprMax.grid.fdtd_grid import FDTDGrid
 from gprMax.grid.mpi_grid import MPIGrid
-from gprMax.output_controllers.geometry_views import GeometryView, GridViewType, Metadata
-from gprMax.output_controllers.grid_view import MPIGridView
+from gprMax.output_controllers.geometry_views import GeometryView, Metadata, MPIMetadata
+from gprMax.output_controllers.grid_view import GridType, MPIGridView
 from gprMax.subgrids.grid import SubGridBaseGrid
 from gprMax.vtkhdf_filehandlers.vtk_unstructured_grid import VtkUnstructuredGrid
 from gprMax.vtkhdf_filehandlers.vtkhdf import VtkCellType
@@ -33,12 +33,27 @@ from gprMax.vtkhdf_filehandlers.vtkhdf import VtkCellType
 logger = logging.getLogger(__name__)
 
 
-class GeometryViewLines(GeometryView[GridViewType]):
+class GeometryViewLines(GeometryView[GridType]):
     """Unstructured grid for a per-cell-edge geometry view."""
 
+    def __init__(
+        self,
+        xs: int,
+        ys: int,
+        zs: int,
+        xf: int,
+        yf: int,
+        zf: int,
+        filename: str,
+        grid: GridType,
+    ):
+        super().__init__(xs, ys, zs, xf, yf, zf, 1, 1, 1, filename, grid)
+
     def prep_vtk(self):
         """Prepares data for writing to VTKHDF file."""
 
+        self.grid_view.initialise_materials(filter_materials=False)
+
         ID = self.grid_view.get_ID()
 
         x = np.arange(self.grid_view.nx + 1, dtype=np.float64)
@@ -81,6 +96,9 @@ class GeometryViewLines(GeometryView[GridViewType]):
             + self.material_data.nbytes
         )
 
+        # Use sorted material IDs rather than default ordering
+        self.material_data = self.grid_view.map_to_view_materials(self.material_data)
+
     def write_vtk(self):
         """Writes geometry information to a VTKHDF file."""
 
@@ -96,21 +114,20 @@ class GeometryViewLines(GeometryView[GridViewType]):
             self.metadata.write_to_vtkhdf(f)
 
 
-class MPIGeometryViewLines(GeometryViewLines[MPIGridView]):
+class MPIGeometryViewLines(GeometryViewLines[MPIGrid]):
     """Image data for a per-cell geometry view."""
 
-    def __init__(self, grid_view: MPIGridView, filename: str, comm: MPI.Comm):
-        super().__init__(grid_view, filename)
-
-        self.comm = comm
-
     @property
-    def grid(self) -> MPIGrid:
-        return self.grid_view.grid
+    def GRID_VIEW_TYPE(self) -> type[MPIGridView]:
+        return MPIGridView
 
     def prep_vtk(self):
         """Prepares data for writing to VTKHDF file."""
 
+        assert isinstance(self.grid_view, MPIGridView)
+
+        self.grid_view.initialise_materials(filter_materials=False)
+
         ID = self.grid_view.get_ID()
 
         x = np.arange(self.grid_view.gx + 1, dtype=np.float64)
@@ -121,11 +138,6 @@ class MPIGeometryViewLines(GeometryViewLines[MPIGridView]):
         self.points += self.grid_view.global_start
         self.points *= self.grid_view.step * self.grid.dl
 
-        # Add offset to subgrid geometry to correctly locate within main grid
-        if isinstance(self.grid, SubGridBaseGrid):
-            offset = [self.grid.i0, self.grid.j0, self.grid.k0]
-            self.points += offset * self.grid.dl * self.grid.ratio
-
         # Each point is the 'source' for 3 lines.
         # NB: Excluding points at the far edge of the geometry as those
         # are the 'source' for no lines
@@ -142,7 +154,7 @@ class MPIGeometryViewLines(GeometryViewLines[MPIGridView]):
         assert isinstance(self.material_data, np.ndarray)
 
         # Write information about any PMLs, sources, receivers
-        self.metadata = Metadata(self.grid_view, averaged_materials=True, materials_only=True)
+        self.metadata = MPIMetadata(self.grid_view, averaged_materials=True, materials_only=True)
 
         # Number of bytes of data to be written to file
         self.nbytes = (
@@ -154,18 +166,20 @@ class MPIGeometryViewLines(GeometryViewLines[MPIGridView]):
         )
 
         # Use global material IDs rather than local IDs
-        self.material_data = self.grid.local_to_global_material_id_map(self.material_data)
+        self.material_data = self.grid_view.map_to_view_materials(self.material_data)
 
     def write_vtk(self):
         """Writes geometry information to a VTKHDF file."""
 
+        assert isinstance(self.grid_view, MPIGridView)
+
         with VtkUnstructuredGrid(
             self.filename,
             self.points,
             self.cell_types,
             self.connectivity,
             self.cell_offsets,
-            comm=self.comm,
+            comm=self.grid_view.comm,
         ) as f:
             f.add_cell_data("Material", self.material_data, self.grid_view.offset)
             self.metadata.write_to_vtkhdf(f)

gprMax/output_controllers/geometry_view_voxels.py

@@ -23,15 +23,15 @@ from mpi4py import MPI
 
 from gprMax._version import __version__
 from gprMax.grid.mpi_grid import MPIGrid
-from gprMax.output_controllers.geometry_views import GeometryView, GridViewType, Metadata
-from gprMax.output_controllers.grid_view import MPIGridView
+from gprMax.output_controllers.geometry_views import GeometryView, Metadata, MPIMetadata
+from gprMax.output_controllers.grid_view import GridType, MPIGridView
 from gprMax.subgrids.grid import SubGridBaseGrid
 from gprMax.vtkhdf_filehandlers.vtk_image_data import VtkImageData
 
 logger = logging.getLogger(__name__)
 
 
-class GeometryViewVoxels(GeometryView[GridViewType]):
+class GeometryViewVoxels(GeometryView[GridType]):
     """Image data for a per-cell geometry view."""
 
     def prep_vtk(self):
@@ -65,31 +65,43 @@ class GeometryViewVoxels(GeometryView[GridViewType]):
             self.metadata.write_to_vtkhdf(f)
 
 
-class MPIGeometryViewVoxels(GeometryViewVoxels[MPIGridView]):
+class MPIGeometryViewVoxels(GeometryViewVoxels[MPIGrid]):
     """Image data for a per-cell geometry view."""
 
-    def __init__(self, grid_view: MPIGridView, filename: str, comm: MPI.Comm):
-        super().__init__(grid_view, filename)
-
-        self.comm = comm
-
     @property
-    def grid(self) -> MPIGrid:
-        return self.grid_view.grid
+    def GRID_VIEW_TYPE(self) -> type[MPIGridView]:
+        return MPIGridView
 
     def prep_vtk(self):
         """Prepares data for writing to VTKHDF file."""
 
-        super().prep_vtk()
+        assert isinstance(self.grid_view, self.GRID_VIEW_TYPE)
+
+        self.material_data = self.grid_view.get_solid()
+
+        self.origin = self.grid_view.global_start * self.grid.dl
+        self.spacing = self.grid_view.step * self.grid.dl
 
         # Use global material IDs rather than local IDs
-        self.material_data = self.grid.local_to_global_material_id_map(self.material_data)
+        self.grid_view.initialise_materials(filter_materials=False)
+        self.material_data = self.grid_view.map_to_view_materials(self.material_data)
+
+        self.nbytes = self.material_data.nbytes
+
+        # Write information about any PMLs, sources, receivers
+        self.metadata = MPIMetadata(self.grid_view)
 
     def write_vtk(self):
         """Writes geometry information to a VTKHDF file."""
 
+        assert isinstance(self.grid_view, self.GRID_VIEW_TYPE)
+
         with VtkImageData(
-            self.filename, self.grid_view.global_size, self.origin, self.spacing, comm=self.comm
+            self.filename,
+            self.grid_view.global_size,
+            self.origin,
+            self.spacing,
+            comm=self.grid_view.comm,
         ) as f:
             f.add_cell_data("Material", self.material_data, self.grid_view.offset)
             self.metadata.write_to_vtkhdf(f)

gprMax/output_controllers/geometry_views.py

@@ -20,17 +20,18 @@ import logging
 import sys
 from abc import abstractmethod
 from pathlib import Path
-from typing import Dict, Generic, List, Sequence, Tuple, Union
+from typing import Dict, Generic, List, Optional, Sequence, Tuple, Union
 
 import h5py
 import numpy as np
 import numpy.typing as npt
+from mpi4py import MPI
 from tqdm import tqdm
 
 import gprMax.config as config
 from gprMax._version import __version__
-from gprMax.grid.fdtd_grid import FDTDGrid
-from gprMax.output_controllers.grid_view import GridType, GridView, GridViewType
+from gprMax.grid.mpi_grid import MPIGrid
+from gprMax.output_controllers.grid_view import GridType, GridView, MPIGridView
 from gprMax.receivers import Rx
 from gprMax.sources import Source
 from gprMax.utilities.utilities import get_terminal_width
@@ -66,15 +67,28 @@ def save_geometry_views(gvs: "List[GeometryView]"):
     logger.info("")
 
 
-class GeometryView(Generic[GridViewType]):
+class GeometryView(Generic[GridType]):
     """Base class for Geometry Views."""
 
     FILE_EXTENSION = ".vtkhdf"
 
+    @property
+    def GRID_VIEW_TYPE(self) -> type[GridView]:
+        return GridView
+
     def __init__(
         self,
-        grid_view: GridViewType,
+        xs: int,
+        ys: int,
+        zs: int,
+        xf: int,
+        yf: int,
+        zf: int,
+        dx: int,
+        dy: int,
+        dz: int,
         filename: str,
+        grid: GridType,
     ):
         """
         Args:
@@ -83,7 +97,7 @@ class GeometryView(Generic[GridViewType]):
             filename: string for filename.
             grid: FDTDGrid class describing a grid in a model.
         """
-        self.grid_view = grid_view
+        self.grid_view = self.GRID_VIEW_TYPE(grid, xs, ys, zs, xf, yf, zf, dx, dy, dz)
 
         self.filenamebase = filename
         self.nbytes = None
@@ -92,7 +106,7 @@ class GeometryView(Generic[GridViewType]):
         self.materials = None
 
     @property
-    def grid(self) -> FDTDGrid:
+    def grid(self) -> GridType:
         return self.grid_view.grid
 
     def set_filename(self):
@@ -118,107 +132,103 @@ class Metadata(Generic[GridType]):
 
     def __init__(
         self,
-        gv: GridView[GridType],
+        grid_view: GridView[GridType],
         averaged_materials: bool = False,
         materials_only: bool = False,
     ):
-        self.gv = gv
+        self.grid_view = grid_view
         self.averaged_materials = averaged_materials
         self.materials_only = materials_only
 
         self.gprmax_version = __version__
-        self.dx_dy_dz = self.grid.dl
-        self.nx_ny_nz = np.array([self.grid.nx, self.grid.ny, self.grid.nz], dtype=np.intc)
+        self.dx_dy_dz = self.dx_dy_dz_comment()
+        self.nx_ny_nz = self.nx_ny_nz_comment()
 
         self.materials = self.materials_comment()
 
         # Write information on PMLs, sources, and receivers
         if not self.materials_only:
             # Information on PML thickness
-            if self.grid.pmls["slabs"]:
-                self.pml_thickness = self.pml_gv_comment()
-            else:
-                self.pml_thickness = None
-            srcs = (
+            self.pml_thickness = self.pml_gv_comment()
+
+            sources = (
                 self.grid.hertziandipoles
                 + self.grid.magneticdipoles
                 + self.grid.voltagesources
                 + self.grid.transmissionlines
             )
-            if srcs:
-                self.source_ids, self.source_positions = self.srcs_rx_gv_comment(srcs)
+            sources_comment = self.srcs_rx_gv_comment(sources)
+            if sources_comment is None:
+                self.source_ids = self.source_positions = None
             else:
-                self.source_ids = None
-                self.source_positions = None
-            if self.grid.rxs:
-                self.receiver_ids, self.receiver_positions = self.srcs_rx_gv_comment(self.grid.rxs)
+                self.source_ids, self.source_positions = sources_comment
+
+            receivers_comment = self.srcs_rx_gv_comment(self.grid.rxs)
+            if receivers_comment is None:
+                self.receiver_ids = self.receiver_positions = None
             else:
-                self.receiver_ids = None
-                self.receiver_positions = None
+                self.receiver_ids, self.receiver_positions = receivers_comment
 
     @property
     def grid(self) -> GridType:
-        return self.gv.grid
+        return self.grid_view.grid
 
     def write_to_vtkhdf(self, file_handler: VtkHdfFile):
-        file_handler.add_field_data(
-            "gprMax_version",
-            self.gprmax_version,
-            dtype=h5py.string_dtype(),
-        )
+        file_handler.add_field_data("gprMax_version", self.gprmax_version)
         file_handler.add_field_data("dx_dy_dz", self.dx_dy_dz)
         file_handler.add_field_data("nx_ny_nz", self.nx_ny_nz)
 
-        file_handler.add_field_data(
-            "material_ids",
-            self.materials,
-            dtype=h5py.string_dtype(),
-        )
+        self.write_material_ids(file_handler)
 
         if not self.materials_only:
             if self.pml_thickness is not None:
                 file_handler.add_field_data("pml_thickness", self.pml_thickness)
 
             if self.source_ids is not None and self.source_positions is not None:
-                file_handler.add_field_data(
-                    "source_ids", self.source_ids, dtype=h5py.string_dtype()
-                )
+                file_handler.add_field_data("source_ids", self.source_ids)
                 file_handler.add_field_data("sources", self.source_positions)
 
             if self.receiver_ids is not None and self.receiver_positions is not None:
-                file_handler.add_field_data(
-                    "receiver_ids", self.receiver_ids, dtype=h5py.string_dtype()
-                )
+                file_handler.add_field_data("receiver_ids", self.receiver_ids)
                 file_handler.add_field_data("receivers", self.receiver_positions)
 
-    def pml_gv_comment(self) -> List[int]:
+    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
 
+        if not grid.pmls["slabs"]:
+            return None
+
         # Only render PMLs if they are in the geometry view
         thickness: Dict[str, int] = grid.pmls["thickness"]
         gv_pml_depth = dict.fromkeys(thickness, 0)
 
-        if self.gv.xs < thickness["x0"]:
-            gv_pml_depth["x0"] = thickness["x0"] - self.gv.xs
-        if self.gv.ys < thickness["y0"]:
-            gv_pml_depth["y0"] = thickness["y0"] - self.gv.ys
-        if thickness["z0"] - self.gv.zs > 0:
-            gv_pml_depth["z0"] = thickness["z0"] - self.gv.zs
-        if self.gv.xf > grid.nx - thickness["xmax"]:
-            gv_pml_depth["xmax"] = self.gv.xf - (grid.nx - thickness["xmax"])
-        if self.gv.yf > grid.ny - thickness["ymax"]:
-            gv_pml_depth["ymax"] = self.gv.yf - (grid.ny - thickness["ymax"])
-        if self.gv.zf > grid.nz - thickness["zmax"]:
-            gv_pml_depth["zmax"] = self.gv.zf - (grid.nz - thickness["zmax"])
+        if self.grid_view.xs < thickness["x0"]:
+            gv_pml_depth["x0"] = thickness["x0"] - self.grid_view.xs
+        if self.grid_view.ys < thickness["y0"]:
+            gv_pml_depth["y0"] = thickness["y0"] - self.grid_view.ys
+        if thickness["z0"] - self.grid_view.zs > 0:
+            gv_pml_depth["z0"] = thickness["z0"] - self.grid_view.zs
+        if self.grid_view.xf > grid.nx - thickness["xmax"]:
+            gv_pml_depth["xmax"] = self.grid_view.xf - (grid.nx - thickness["xmax"])
+        if self.grid_view.yf > grid.ny - thickness["ymax"]:
+            gv_pml_depth["ymax"] = self.grid_view.yf - (grid.ny - thickness["ymax"])
+        if self.grid_view.zf > grid.nz - thickness["zmax"]:
+            gv_pml_depth["zmax"] = self.grid_view.zf - (grid.nz - thickness["zmax"])
 
-        return list(gv_pml_depth.values())
+        return np.array(list(gv_pml_depth.values()), dtype=np.int64)
 
     def srcs_rx_gv_comment(
         self, srcs: Union[Sequence[Source], List[Rx]]
-    ) -> Tuple[List[str], npt.NDArray[np.float32]]:
+    ) -> Optional[Tuple[List[str], npt.NDArray[np.float64]]]:
         """Used to name sources and/or receivers."""
+        if not srcs:
+            return None
+
         names: List[str] = []
-        positions: npt.NDArray[np.float32] = np.empty((len(srcs), 3))
+        positions = np.empty((len(srcs), 3))
         for index, src in enumerate(srcs):
             position = src.coord * self.grid.dl
             names.append(src.ID)
@@ -229,11 +239,69 @@ class Metadata(Generic[GridType]):
     def dx_dy_dz_comment(self) -> npt.NDArray[np.float64]:
         return self.grid.dl
 
-    def nx_ny_nz_comment(self) -> npt.NDArray[np.intc]:
-        return np.array([self.grid.nx, self.grid.ny, self.grid.nz], dtype=np.intc)
+    def nx_ny_nz_comment(self) -> npt.NDArray[np.int32]:
+        return self.grid.size
+
+    def materials_comment(self) -> Optional[List[str]]:
+        if self.grid_view.materials is None:
+            return None
 
-    def materials_comment(self) -> List[str]:
         if not self.averaged_materials:
-            return [m.ID for m in self.grid.materials if m.type != "dielectric-smoothed"]
+            return [m.ID for m in self.grid_view.materials if m.type != "dielectric-smoothed"]
         else:
-            return [m.ID for m in self.grid.materials]
+            return [m.ID for m in self.grid_view.materials]
+
+
+class MPIMetadata(Metadata[MPIGrid]):
+    def nx_ny_nz_comment(self) -> npt.NDArray[np.int32]:
+        return self.grid.global_size
+
+    def pml_gv_comment(self) -> Optional[npt.NDArray[np.int64]]:
+        gv_pml_depth = super().pml_gv_comment()
+
+        if gv_pml_depth is None:
+            gv_pml_depth = np.zeros(6, dtype=np.int64)
+
+        assert isinstance(self.grid_view, MPIGridView)
+        recv_buffer = np.empty((self.grid_view.comm.size, 6), dtype=np.int64)
+        self.grid_view.comm.Allgather(gv_pml_depth, recv_buffer)
+
+        gv_pml_depth = np.max(recv_buffer, axis=0)
+
+        return None if all(gv_pml_depth == 0) else gv_pml_depth
+
+    def srcs_rx_gv_comment(
+        self, srcs: Union[Sequence[Source], List[Rx]]
+    ) -> Optional[Tuple[List[str], npt.NDArray[np.float64]]]:
+        objects: Dict[str, npt.NDArray[np.float64]] = {}
+        for src in srcs:
+            position = self.grid.local_to_global_coordinate(src.coord) * self.grid.dl
+            objects[src.ID] = position
+
+        assert isinstance(self.grid_view, MPIGridView)
+        global_objects: List[Dict[str, npt.NDArray[np.float64]]] = self.grid_view.comm.allgather(
+            objects
+        )
+        objects = {k: v for d in global_objects for k, v in d.items()}
+        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)

gprMax/output_controllers/grid_view.py

@@ -1,6 +1,6 @@
 import logging
 from itertools import chain
-from typing import Generic
+from typing import Generic, Tuple
 
 import numpy as np
 import numpy.typing as npt
@@ -62,6 +62,42 @@ class GridView(Generic[GridType]):
             f" step={self.step}, size={self.size})"
         )
 
+    @property
+    def xs(self) -> int:
+        return self.start[0]
+
+    @property
+    def ys(self) -> int:
+        return self.start[1]
+
+    @property
+    def zs(self) -> int:
+        return self.start[2]
+
+    @property
+    def xf(self) -> int:
+        return self.stop[0]
+
+    @property
+    def yf(self) -> int:
+        return self.stop[1]
+
+    @property
+    def zf(self) -> int:
+        return self.stop[2]
+
+    @property
+    def nx(self) -> int:
+        return self.size[0]
+
+    @property
+    def ny(self) -> int:
+        return self.size[1]
+
+    @property
+    def nz(self) -> int:
+        return self.size[2]
+
     def get_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
         """Create a slice object for the specified dimension.
 
@@ -112,7 +148,7 @@ class GridView(Generic[GridType]):
             ]
         )
 
-    def initialise_materials(self):
+    def initialise_materials(self, filter_materials: bool = True):
         """Create a new ID map for materials in the grid view.
 
         Rather than using the default material IDs (as per the main grid
@@ -125,11 +161,16 @@ class GridView(Generic[GridType]):
         map_to_view_materials() function.
         """
         # Get unique materials in the grid view
-        ID = self.get_ID(force_refresh=True)
-        materials_in_grid_view = np.unique(ID)
+        if filter_materials:
+            ID = self.get_ID(force_refresh=True)
+            materials_in_grid_view = np.unique(ID)
 
-        # Get actual Material objects and sort
-        self.materials = np.array(self.grid.materials, dtype=Material)[materials_in_grid_view]
+            # Get actual Material objects
+            self.materials = np.array(self.grid.materials, dtype=Material)[materials_in_grid_view]
+        else:
+            self.materials = np.array(self.grid.materials, dtype=Material)
+
+        # Sort materials
         self.materials.sort()
 
         # Create map from material ID to 0 - number of materials
@@ -150,7 +191,7 @@ class GridView(Generic[GridType]):
         Returns:
             array: Mapped array.
         """
-        return self.map_materials_func(array)
+        return self.map_materials_func(array).astype(array.dtype)
 
     def get_ID(self, force_refresh=False) -> npt.NDArray[np.uint32]:
         """Get a view of the ID array.
@@ -277,6 +318,12 @@ class MPIGridView(GridView[MPIGrid]):
         """
         super().__init__(grid, xs, ys, zs, xf, yf, zf, dx, dy, dz)
 
+        comm = grid.comm.Split()
+        assert isinstance(comm, MPI.Intracomm)
+        start_grid_coord = grid.get_grid_coord_from_local_coordinate(self.start)
+        stop_grid_coord = grid.get_grid_coord_from_local_coordinate(self.stop) + 1
+        self.comm = comm.Create_cart((stop_grid_coord - start_grid_coord).tolist())
+
         self.global_size = self.size
 
         # Calculate start for the local grid
@@ -324,6 +371,18 @@ class MPIGridView(GridView[MPIGrid]):
             f" stop={self.stop}, step={self.step}, size={self.size}, offset={self.offset})"
         )
 
+    @property
+    def gx(self) -> int:
+        return self.global_size[0]
+
+    @property
+    def gy(self) -> int:
+        return self.global_size[1]
+
+    @property
+    def gz(self) -> int:
+        return self.global_size[2]
+
     def get_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
         """Create a slice object for the specified dimension.
 
@@ -348,6 +407,35 @@ class MPIGridView(GridView[MPIGrid]):
         return slice(self.start[dimension], stop, self.step[dimension])
 
     def get_output_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
+        """Create an output slice object for the specified dimension.
+
+        This provides a slice of the grid view for the section of the
+        grid view managed by this rank. This can be used when writing
+        out arrays provided by the grid view as part of a collective
+        operation.
+
+        For example:
+        ```
+        dset_slice = (
+                grid_view.get_output_slice(0),
+                grid_view.get_output_slice(1),
+                grid_view.get_output_slice(2),
+            )
+
+        dset[dset_slice] = grid_view.get_solid()
+        ```
+
+        Args:
+            dimension: Dimension to create the slice object for. Values
+                0, 1, and 2 map to the x, y, and z dimensions
+                respectively.
+            upper_bound_exclusive: Optionally specify if the upper bound
+                of the slice should be exclusive or inclusive. Defaults
+                to True.
+
+        Returns:
+            slice: Slice object
+        """
         if upper_bound_exclusive or self.has_positive_neighbour[dimension]:
             size = self.size[dimension]
         else:
@@ -359,7 +447,32 @@ class MPIGridView(GridView[MPIGrid]):
 
         return slice(offset, offset + size)
 
-    def initialise_materials(self, comm: MPI.Cartcomm):
+    def get_3d_output_slice(self, upper_bound_exclusive: bool = True) -> Tuple[slice, slice, slice]:
+        """Create a 3D output slice object.
+
+        This provides a slice of the grid view for the section of the
+        grid view managed by this rank. This can be used when writing
+        out arrays provided by the grid view as part of a collective
+        operation.
+
+        For example:
+        `dset[grid_view.get_3d_output_slice()] = grid_view.get_solid()`
+
+        Args:
+            upper_bound_exclusive: Optionally specify if the upper bound
+                of the slice should be exclusive or inclusive. Defaults
+                to True.
+
+        Returns:
+            slice: 3D Slice object
+        """
+        return (
+            self.get_output_slice(0, upper_bound_exclusive),
+            self.get_output_slice(1, upper_bound_exclusive),
+            self.get_output_slice(2, upper_bound_exclusive),
+        )
+
+    def initialise_materials(self, filter_materials: bool = True):
         """Create a new ID map for materials in the grid view.
 
         Rather than using the default material IDs (as per the main grid
@@ -372,15 +485,18 @@ class MPIGridView(GridView[MPIGrid]):
         communication to construct the new map. It should also be called
         before the map_to_view_materials() function.
         """
-        ID = self.get_ID(force_refresh=True)
+        if filter_materials:
+            ID = self.get_ID(force_refresh=True)
+            local_material_ids = np.unique(ID)
+            local_materials = np.array(self.grid.materials, dtype=Material)[local_material_ids]
+        else:
+            local_materials = np.array(self.grid.materials, dtype=Material)
 
-        local_material_ids = np.unique(ID)
-        local_materials = np.array(self.grid.materials, dtype=Material)[local_material_ids]
         local_materials.sort()
         local_material_ids = [m.numID for m in local_materials]
 
         # Send all materials to the coordinating rank
-        materials_by_rank = comm.gather(local_materials, root=0)
+        materials_by_rank = self.comm.gather(local_materials, root=0)
 
         if materials_by_rank is not None:
             # Filter out duplicate materials and sort by material ID
@@ -390,9 +506,9 @@ class MPIGridView(GridView[MPIGrid]):
             # The new material IDs corespond to each material's index in
             # the sorted self.materials array. For each rank, get the
             # new IDs of each material it sent to send back
-            for rank in range(1, comm.size):
+            for rank in range(1, self.comm.size):
                 new_material_ids = np.where(np.isin(self.materials, materials_by_rank[rank]))[0]
-                comm.Isend([new_material_ids.astype(np.int32), MPI.INT], rank)
+                self.comm.Isend([new_material_ids.astype(np.int32), MPI.INT], rank)
 
             new_material_ids = np.where(np.isin(self.materials, materials_by_rank[0]))[0]
             new_material_ids = new_material_ids.astype(np.int32)
@@ -401,7 +517,7 @@ class MPIGridView(GridView[MPIGrid]):
 
             # Get list of global IDs for this rank's local materials
             new_material_ids = np.empty(len(local_materials), dtype=np.int32)
-            comm.Recv([new_material_ids, MPI.INT], 0)
+            self.comm.Recv([new_material_ids, MPI.INT], 0)
 
         # Create map from local material ID to global material ID
         materials_map = {

gprMax/user_objects/cmds_multiuse.py

@@ -389,16 +389,9 @@ class VoltageSource(RotatableMixin, GridUserObject):
         voltage_source = VoltageSourceUser()
         voltage_source.polarisation = self.polarisation
         voltage_source.coord = coord
-        voltage_source.ID = (
-            voltage_source.__class__.__name__
-            + "("
-            + str(voltage_source.xcoord)
-            + ","
-            + str(voltage_source.ycoord)
-            + ","
-            + str(voltage_source.zcoord)
-            + ")"
-        )
+        uip = self._create_uip(grid)
+        x, y, z = uip.discretise_static_point(self.point)
+        voltage_source.ID = f"{voltage_source.__class__.__name__}({x},{y},{z})"
         voltage_source.resistance = self.resistance
         voltage_source.waveform = grid.get_waveform_by_id(self.waveform_id)
 
@@ -540,7 +533,9 @@ class HertzianDipole(RotatableMixin, GridUserObject):
 
         h.coord = coord
         h.coordorigin = coord
-        h.ID = f"{h.__class__.__name__}({str(h.xcoord)},{str(h.ycoord)},{str(h.zcoord)})"
+        uip = self._create_uip(grid)
+        x, y, z = uip.discretise_static_point(self.point)
+        h.ID = f"{h.__class__.__name__}({x},{y},{z})"
         h.waveform = grid.get_waveform_by_id(self.waveform_id)
 
         if self.start is None or self.stop is None:
@@ -696,16 +691,9 @@ class MagneticDipole(RotatableMixin, GridUserObject):
         m.polarisation = self.polarisation
         m.coord = coord
         m.coordorigin = coord
-        m.ID = (
-            m.__class__.__name__
-            + "("
-            + str(m.xcoord)
-            + ","
-            + str(m.ycoord)
-            + ","
-            + str(m.zcoord)
-            + ")"
-        )
+        uip = self._create_uip(grid)
+        x, y, z = uip.discretise_static_point(self.point)
+        m.ID = f"{m.__class__.__name__}({x},{y},{z})"
         m.waveform = grid.get_waveform_by_id(self.waveform_id)
 
         if self.start is None or self.stop is None:
@@ -858,16 +846,9 @@ class TransmissionLine(RotatableMixin, GridUserObject):
         t = TransmissionLineUser(grid.iterations, grid.dt)
         t.polarisation = self.polarisation
         t.coord = coord
-        t.ID = (
-            t.__class__.__name__
-            + "("
-            + str(t.xcoord)
-            + ","
-            + str(t.ycoord)
-            + ","
-            + str(t.zcoord)
-            + ")"
-        )
+        uip = self._create_uip(grid)
+        x, y, z = uip.discretise_static_point(self.point)
+        t.ID = f"{t.__class__.__name__}({x},{y},{z})"
         t.resistance = self.resistance
         t.waveform = grid.get_waveform_by_id(self.waveform_id)
 

gprMax/user_objects/cmds_output.py

@@ -34,8 +34,20 @@ class GeometryView(OutputUserObject):
     def hash(self):
         return "#geometry_view"
 
-    def __init__(self, **kwargs):
-        super().__init__(**kwargs)
+    def __init__(
+        self,
+        p1: Tuple[float, float, float],
+        p2: Tuple[float, float, float],
+        dl: Tuple[float, float, float],
+        output_type: str,
+        filename: str,
+    ):
+        super().__init__(p1=p1, p2=p2, dl=dl, filename=filename, output_type=output_type)
+        self.lower_bound = p1
+        self.upper_bound = p2
+        self.dl = dl
+        self.filename = filename
+        self.output_type = output_type
 
     def geometry_view_constructor(self, output_type):
         """Selects appropriate class for geometry view dependent on geometry
@@ -50,71 +62,61 @@ class GeometryView(OutputUserObject):
         return GeometryViewUser
 
     def build(self, model: Model, grid: FDTDGrid):
-        try:
-            p1 = self.kwargs["p1"]
-            p2 = self.kwargs["p2"]
-            dl = self.kwargs["dl"]
-            output_type = self.kwargs["output_type"].lower()
-            filename = self.kwargs["filename"]
-        except KeyError:
-            logger.exception(f"{self.params_str()} requires exactly eleven parameters.")
-            raise
-
-        GeometryViewUser = self.geometry_view_constructor(output_type)
-
         uip = self._create_uip(grid)
-        try:
-            p3 = uip.round_to_grid_static_point(p1)
-            p4 = uip.round_to_grid_static_point(p2)
-            p1, p2 = uip.check_box_points(p1, p2, self.params_str())
-        except ValueError:
-            logger.exception(f"{self.params_str()} point is outside the domain.")
-            raise
-        xs, ys, zs = p1
-        xf, yf, zf = p2
+        discretised_lower_bound, discretised_upper_bound = uip.check_output_object_bounds(
+            self.lower_bound, self.upper_bound, self.params_str()
+        )
+        discretised_dl = uip.discretise_static_point(self.dl)
 
-        dx, dy, dz = uip.discretise_static_point(dl)
-
-        if dx < 0 or dy < 0 or dz < 0:
-            logger.exception(f"{self.params_str()} the step size should not be less than zero.")
-            raise ValueError
-        if dx > grid.nx or dy > grid.ny or dz > grid.nz:
-            logger.exception(
+        if any(discretised_dl < 0):
+            raise ValueError(f"{self.params_str()} the step size should not be less than zero.")
+        if any(discretised_dl > grid.size):
+            raise ValueError(
                 f"{self.params_str()} the step size should be less than the domain size."
             )
-            raise ValueError
-        if dx < 1 or dy < 1 or dz < 1:
-            logger.exception(
-                f"{self.params_str()} the step size should not be less than the spatial discretisation."
+        if any(discretised_dl < 1):
+            raise ValueError(
+                f"{self.params_str()} the step size should not be less than the spatial"
+                " discretisation."
             )
-            raise ValueError
-        if output_type not in ["n", "f"]:
-            logger.exception(
+        if self.output_type == "f" and any(discretised_dl != 1):
+            raise ValueError(
+                f"{self.params_str()} requires the spatial discretisation for the geometry view to"
+                " be the same as the model for geometry view of type f (fine)."
+            )
+
+        if self.output_type == "n":
+            g = model.add_geometry_view_voxels(
+                grid,
+                discretised_lower_bound,
+                discretised_upper_bound,
+                discretised_dl,
+                self.filename,
+            )
+        elif self.output_type == "f":
+            g = model.add_geometry_view_lines(
+                grid,
+                discretised_lower_bound,
+                discretised_upper_bound,
+                self.filename,
+            )
+        else:
+            raise ValueError(
                 f"{self.params_str()} requires type to be either n (normal) or f (fine)."
             )
-            raise ValueError
-        if output_type == "f" and (
-            dx * grid.dx != grid.dx or dy * grid.dy != grid.dy or dz * grid.dz != grid.dz
-        ):
-            logger.exception(
-                f"{self.params_str()} requires the spatial "
-                "discretisation for the geometry view to be the "
-                "same as the model for geometry view of "
-                "type f (fine)"
+
+        if g is not None:
+            p1 = uip.round_to_grid_static_point(self.lower_bound)
+            p2 = uip.round_to_grid_static_point(self.upper_bound)
+            dl = discretised_dl * grid.dl
+
+            logger.info(
+                f"{self.grid_name(grid)}Geometry view from"
+                f" {p1[0]:g}m, {p1[1]:g}m, {p1[2]:g}m,"
+                f" to {p2[0]:g}m, {p2[1]:g}m, {p2[2]:g}m,"
+                f" discretisation {dl[0]:g}m, {dl[1]:g}m, {dl[2]:g}m,"
+                f" with filename base {g.filenamebase} created."
             )
-            raise ValueError
-
-        g = GeometryViewUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, grid)
-
-        logger.info(
-            f"{self.grid_name(grid)}Geometry view from {p3[0]:g}m, "
-            f"{p3[1]:g}m, {p3[2]:g}m, to {p4[0]:g}m, {p4[1]:g}m, "
-            f"{p4[2]:g}m, discretisation {dx * grid.dx:g}m, "
-            f"{dy * grid.dy:g}m, {dz * grid.dz:g}m, with filename "
-            f"base {g.filename} created."
-        )
-
-        model.geometryviews.append(g)
 
 
 class GeometryObjectsWrite(OutputUserObject):

gprMax/vtkhdf_filehandlers/vtkhdf.py

@@ -5,7 +5,7 @@ from enum import Enum
 from os import PathLike
 from pathlib import Path
 from types import TracebackType
-from typing import Optional, Union
+from typing import Optional, Tuple, Union
 
 import h5py
 import numpy as np
@@ -211,7 +211,7 @@ class VtkHdfFile(AbstractContextManager):
         self,
         path: str,
         data: npt.ArrayLike,
-        shape: Optional[npt.NDArray[np.int32]] = None,
+        shape: Optional[Union[npt.NDArray[np.int32], Tuple[int, ...]]] = None,
         offset: Optional[npt.NDArray[np.int32]] = None,
         dtype: Optional[npt.DTypeLike] = None,
         xyz_data_ordering=True,
@@ -228,6 +228,8 @@ class VtkHdfFile(AbstractContextManager):
                 Can be omitted if data provides the full dataset.
             offset (optional): Offset to store the provided data at. Can
                 be omitted if data provides the full dataset.
+            dtype (optional): Type of the data. If omitted, the type
+                will be deduced from the provided data.
             xyz_data_ordering (optional): If True, the data will be
                 transposed as VTKHDF stores datasets using ZYX ordering.
                 Default True.
@@ -237,17 +239,37 @@ class VtkHdfFile(AbstractContextManager):
                 and offset are invalid.
         """
 
-        # If dtype is a string, ensure it is ascii encoded and use
-        # variable length strings
-        if dtype is not None and h5py.check_string_dtype(dtype) is not None:
-            dtype = h5py.string_dtype(encoding="ascii", length=None)
+        # If dtype is a string and using parallel I/O, ensure using
+        # fixed length strings
+        if isinstance(dtype, np.dtype) and self.comm is not None:
+            string_info = h5py.check_string_dtype(dtype)
+            if string_info is not None and string_info.length is None:
+                logger.warning(
+                    "HDF5 does not support variable length strings with parallel I/O."
+                    " Using fixed length strings instead."
+                )
+                dtype = h5py.string_dtype(encoding="ascii", length=0)
 
         if not isinstance(data, np.ndarray):
             data = np.array(data, dtype=dtype)
             if data.ndim < 1:
                 data = np.expand_dims(data, axis=-1)
 
-        if dtype is None:
+        if data.dtype.kind == "U":
+            if dtype is not None:  # Only log warning if user specified a data type
+                logger.warning(
+                    "NumPy UTF-32 ('U' dtype) is not supported by HDF5."
+                    " Converting to bytes array ('S' dtype)."
+                )
+            data = data.astype("S")
+
+        # Explicitly define string datatype
+        # VTKHDF only supports ascii strings (not UTF-8)
+        if data.dtype.kind == "S":
+            dtype = h5py.string_dtype(encoding="ascii", length=data.dtype.itemsize)
+            data = data.astype(dtype)
+
+        elif dtype is None:
             dtype = data.dtype
 
         # VTKHDF stores datasets using ZYX ordering rather than XYZ
@@ -260,8 +282,14 @@ class VtkHdfFile(AbstractContextManager):
         if offset is not None:
             offset = np.flip(offset)
 
+        logger.debug(
+            f"Writing dataset '{path}', shape: {shape}, data.shape: {data.shape}, dtype: {dtype}"
+        )
+
         if shape is None or all(shape == data.shape):
-            self.file_handler.create_dataset(path, data=data, dtype=dtype)
+            shape = data.shape if shape is None else shape
+            dataset = self.file_handler.create_dataset(path, shape=shape, dtype=dtype)
+            dataset[:] = data
         elif offset is None:
             raise ValueError(
                 "Offset must not be None as the full dataset has not been provided."
@@ -293,15 +321,56 @@ class VtkHdfFile(AbstractContextManager):
             start = offset
             stop = offset + data.shape
 
-            dataset_slice = (slice(start[i], stop[i]) for i in range(dimensions))
+            dataset_slice = tuple([slice(start[i], stop[i]) for i in range(dimensions)])
 
             dataset[dataset_slice] = data
 
+    def _create_dataset(
+        self, path: str, shape: Union[npt.NDArray[np.int32], Tuple[int, ...]], dtype: npt.DTypeLike
+    ):
+        """Create dataset in the VTKHDF file without writing any data.
+
+        Args:
+            path: Absolute path to the dataset.
+            shape: Size of the full dataset being created.
+            dtype: Type of the data.
+
+        Raises:
+            TypeError: Raised if attempt to use variable length strings
+                with parallel I/O.
+        """
+        dtype = np.dtype(dtype)
+
+        # If dtype is a string and using parallel I/O, ensure using
+        # fixed length strings
+        if self.comm is not None:
+            string_info = h5py.check_string_dtype(dtype)
+            if string_info is not None and string_info.length is None:
+                raise TypeError(
+                    "HDF5 does not support variable length strings with parallel I/O."
+                    " Use fixed length strings instead."
+                )
+
+        if dtype.kind == "U":
+            logger.warning(
+                "NumPy UTF-32 ('U' dtype) is not supported by HDF5."
+                " Converting to bytes array ('S' dtype)."
+            )
+
+        # Explicitly define string datatype
+        # VTKHDF only supports ascii strings (not UTF-8)
+        if dtype.kind == "U" or dtype.kind == "S":
+            dtype = h5py.string_dtype(encoding="ascii", length=dtype.itemsize)
+
+        logger.debug(f"Creating dataset '{path}', shape: {shape}, dtype: {dtype}")
+
+        self.file_handler.create_dataset(path, shape=shape, dtype=dtype)
+
     def add_point_data(
         self,
         name: str,
         data: npt.NDArray,
-        shape: Optional[npt.NDArray[np.int32]] = None,
+        shape: Optional[Union[npt.NDArray[np.int32], Tuple[int, ...]]] = None,
         offset: Optional[npt.NDArray[np.int32]] = None,
     ):
         """Add point data to the VTKHDF file.
@@ -321,7 +390,7 @@ class VtkHdfFile(AbstractContextManager):
         self,
         name: str,
         data: npt.NDArray,
-        shape: Optional[npt.NDArray[np.int32]] = None,
+        shape: Optional[Union[npt.NDArray[np.int32], Tuple[int, ...]]] = None,
         offset: Optional[npt.NDArray[np.int32]] = None,
     ):
         """Add cell data to the VTKHDF file.
@@ -340,8 +409,8 @@ class VtkHdfFile(AbstractContextManager):
     def add_field_data(
         self,
         name: str,
-        data: npt.ArrayLike,
-        shape: Optional[npt.NDArray[np.int32]] = None,
+        data: Optional[npt.ArrayLike],
+        shape: Optional[Union[npt.NDArray[np.int32], Tuple[int, ...]]] = None,
         offset: Optional[npt.NDArray[np.int32]] = None,
         dtype: Optional[npt.DTypeLike] = None,
     ):
@@ -349,16 +418,32 @@ class VtkHdfFile(AbstractContextManager):
 
         Args:
             name: Name of the dataset.
-            data: Data to be saved.
+            data: Data to be saved. Can be None if both shape and dtype
+                are specified. If None, the dataset will be created but
+                no data written. This can be useful if, for example,
+                only one rank is writing the data. As long as all ranks
+                know the shape and dtype, ranks not writing data can
+                perform the collective operation of creating the
+                dataset, but only the rank(s) with the data need to
+                write data.
             shape (optional): Size of the full dataset being created.
                 Can be omitted if data provides the full dataset.
             offset (optional): Offset to store the provided data at. Can
                 be omitted if data provides the full dataset.
+            dtype (optional): Type of the data. If omitted, the type
+                will be deduced from the provided data.
         """
         dataset_path = self._build_dataset_path("FieldData", name)
-        self._write_dataset(
-            dataset_path, data, shape=shape, offset=offset, dtype=dtype, xyz_data_ordering=False
-        )
+        if data is not None:
+            self._write_dataset(
+                dataset_path, data, shape=shape, offset=offset, dtype=dtype, xyz_data_ordering=False
+            )
+        elif shape is not None and dtype is not None:
+            self._create_dataset(dataset_path, shape, dtype)
+        else:
+            raise ValueError(
+                "If data is None, shape and dtype must be provided. I.e. they must not be None"
+            )
 
 
 class VtkCellType(np.uint8, Enum):