Fix GeometryObjects for parallel build

Update GeometryObject class to use new GridView class.

gprMax/grid/mpi_grid.py

@@ -138,6 +138,20 @@ class MPIGrid(FDTDGrid):
         """
         return self.rank == self.COORDINATOR_RANK
 
+    def create_sub_communicator(
+        self, start: npt.NDArray[np.int32], stop: npt.NDArray[np.int32]
+    ) -> Optional[MPI.Cartcomm]:
+        if self.local_bounds_overlap_grid(start, stop):
+            comm = self.comm.Split()
+            assert isinstance(comm, MPI.Intracomm)
+            start_grid_coord = self.get_grid_coord_from_coordinate(start)
+            stop_grid_coord = self.get_grid_coord_from_coordinate(stop) + 1
+            comm = comm.Create_cart((stop_grid_coord - start_grid_coord).tolist())
+            return comm
+        else:
+            self.comm.Split(MPI.UNDEFINED)
+            return None
+
     def get_grid_coord_from_coordinate(self, coord: npt.NDArray[np.int32]) -> npt.NDArray[np.int32]:
         """Get the MPI grid coordinate for a global grid coordinate.
 
@@ -252,11 +266,9 @@ class MPIGrid(FDTDGrid):
 
         return all(global_coord >= lower_bound) and all(global_coord <= upper_bound)
 
-    def global_bounds_overlap_local_grid(
-        self, start: npt.NDArray[np.int32], stop: npt.NDArray[np.int32]
+    def local_bounds_overlap_grid(
+        self, local_start: npt.NDArray[np.int32], local_stop: npt.NDArray[np.int32]
     ) -> bool:
-        local_start = self.global_to_local_coordinate(start)
-        local_stop = self.global_to_local_coordinate(stop)
         return all(local_start < self.size) and all(local_stop > self.negative_halo_offset)
 
     def limit_global_bounds_to_within_local_grid(

gprMax/model.py

@@ -19,7 +19,7 @@
 import datetime
 import logging
 import sys
-from typing import List, Sequence, Tuple
+from typing import List, Optional, Sequence
 
 import humanize
 import numpy as np
@@ -29,6 +29,7 @@ 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.subgrids.grid import SubGridBaseGrid
 
 init()
@@ -38,7 +39,7 @@ from tqdm import tqdm
 import gprMax.config as config
 
 from .fields_outputs import write_hdf5_outputfile
-from .geometry_outputs import GeometryObjects, GeometryView, save_geometry_views
+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
@@ -64,7 +65,7 @@ class Model:
         self.subgrids: List[SubGridBaseGrid] = []
 
         self.geometryviews: List[GeometryView] = []
-        self.geometryobjects: List[GeometryObjects] = []
+        self.geometryobjects: List[GeometryObject] = []
 
         # Monitor memory usage
         self.p = None
@@ -173,6 +174,19 @@ class Model:
     def set_size(self, size: npt.NDArray[np.int32]):
         self.nx, self.ny, self.nz = size
 
+    def add_geometry_object(
+        self,
+        grid: FDTDGrid,
+        start: npt.NDArray[np.int32],
+        stop: npt.NDArray[np.int32],
+        basefilename: str,
+    ) -> Optional[GeometryObject]:
+        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 build(self):
         """Builds the Yee cells for a model."""
 
@@ -222,7 +236,7 @@ class Model:
                     file=sys.stdout,
                     disable=not config.sim_config.general["progressbars"],
                 )
-                go.write_hdf5(self.title, self.G, pbar)
+                go.write_hdf5(self.title, pbar)
                 pbar.close()
             logger.info("")
 

gprMax/mpi_model.py

@@ -9,6 +9,7 @@ from gprMax import config
 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.snapshots import save_snapshots
 
 logger = logging.getLogger(__name__)
@@ -59,34 +60,23 @@ class MPIModel(Model):
 
         self.G.calculate_local_extents()
 
-    def build_geometry(self):
-        super().build_geometry()
+    def add_geometry_object(
+        self,
+        grid: MPIGrid,
+        start: npt.NDArray[np.int32],
+        stop: npt.NDArray[np.int32],
+        basefilename: str,
+    ) -> Optional[MPIGeometryObject]:
+        comm = grid.create_sub_communicator(start, stop)
 
-        self._filter_geometry_objects()
-
-    def _filter_geometry_objects(self):
-        objects = self.comm.bcast(self.geometryobjects)
-        self.geometryobjects = []
-
-        for go in objects:
-            start = np.array([go.xs, go.ys, go.zs], dtype=np.int32)
-            stop = np.array([go.xf, go.yf, go.zf], dtype=np.int32)
-            if self.G.global_bounds_overlap_local_grid(start, stop):
-                comm = self.comm.Split()
-                assert isinstance(comm, MPI.Intracomm)
-                start_grid_coord = self.G.get_grid_coord_from_coordinate(start)
-                stop_grid_coord = self.G.get_grid_coord_from_coordinate(stop) + 1
-                go.comm = comm.Create_cart((stop_grid_coord - start_grid_coord).tolist())
-
-                go.global_size = np.array([go.nx, go.ny, go.nz], dtype=np.int32)
-                start, stop, offset = self.G.limit_global_bounds_to_within_local_grid(start, stop)
-                go.size = stop - start
-                go.start = start
-                go.stop = stop
-                go.offset = offset
-                self.geometryobjects.append(go)
-            else:
-                self.comm.Split(MPI.UNDEFINED)
+        if comm is None:
+            return None
+        else:
+            geometry_object = MPIGeometryObject(
+                grid, start[0], start[1], start[2], stop[0], stop[1], stop[2], basefilename, comm
+            )
+            self.geometryobjects.append(geometry_object)
+            return geometry_object
 
     def write_output_data(self):
         """Writes output data, i.e. field data for receivers and snapshots to

gprMax/output_controllers/geometry_objects.py

@@ -10,32 +10,39 @@ from tqdm import tqdm
 from gprMax import config
 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 GridViewType, MPIGridView
+from gprMax.output_controllers.grid_view import GridView, MPIGridView
 
 
-class GeometryObjects(Generic[GridViewType]):
+class GeometryObject:
     """Geometry objects to be written to file."""
 
-    def __init__(self, grid_view: GridViewType, filename: str):
+    @property
+    def GRID_VIEW_TYPE(self) -> type[GridView]:
+        return GridView
+
+    def __init__(
+        self, grid: FDTDGrid, xs: int, ys: int, zs: int, xf: int, yf: int, zf: int, filename: str
+    ):
         """
         Args:
             xs, xf, ys, yf, zs, zf: ints for extent of the volume in cells.
             filename: string for filename.
         """
-        self.grid_view = grid_view
+        self.grid_view = self.GRID_VIEW_TYPE(grid, xs, ys, zs, xf, yf, zf)
 
         # Set filenames
         parts = config.sim_config.input_file_path.with_suffix("").parts
         self.filename = Path(*parts[:-1], filename)
         self.filename_hdf5 = self.filename.with_suffix(".h5")
-        self.filename_materials = Path(self.filename, f"{self.filename}_materials")
+        self.filename_materials = Path(f"{self.filename}_materials")
         self.filename_materials = self.filename_materials.with_suffix(".txt")
 
         # Sizes of arrays to write necessary to update progress bar
-        self.solidsize = np.prod(self.grid_view.size + 1) * np.dtype(np.uint32).itemsize
-        self.rigidsize = 18 * np.prod(self.grid_view.size + 1) * np.dtype(np.int8).itemsize
-        self.IDsize = 6 * np.prod(self.grid_view.size + 1) * np.dtype(np.uint32).itemsize
+        self.solidsize = (float)(np.prod(self.grid_view.size + 1) * np.dtype(np.uint32).itemsize)
+        self.rigidsize = (float)(18 * np.prod(self.grid_view.size + 1) * np.dtype(np.int8).itemsize)
+        self.IDsize = (float)(6 * np.prod(self.grid_view.size + 1) * np.dtype(np.uint32).itemsize)
         self.datawritesize = self.solidsize + self.rigidsize + self.IDsize
 
     @property
@@ -72,7 +79,7 @@ class GeometryObjects(Generic[GridViewType]):
             dispersionstr += material.ID
             file.write(dispersionstr + "\n")
 
-    def write_hdf5(self, title: str, G: FDTDGrid, pbar: tqdm):
+    def write_hdf5(self, title: str, pbar: tqdm):
         """Writes a geometry objects file in HDF5 format.
 
         Args:
@@ -109,15 +116,24 @@ class GeometryObjects(Generic[GridViewType]):
                 self.output_material(material, fmaterials)
 
 
-class MPIGeometryObjects(GeometryObjects[MPIGridView]):
+class MPIGeometryObject(GeometryObject):
+    @property
+    def GRID_VIEW_TYPE(self) -> type[MPIGridView]:
+        return MPIGridView
+
     def __init__(
         self,
-        grid_view: MPIGridView,
+        grid: MPIGrid,
+        xs: int,
+        xf: int,
+        ys: int,
+        yf: int,
+        zs: int,
+        zf: int,
         filename: str,
-        comm: MPI.Comm = None,
+        comm: MPI.Cartcomm,
     ):
-        super().__init__(grid_view, filename)
-
+        super().__init__(grid, xs, xf, ys, yf, zs, zf, filename)
         self.comm = comm
 
     def write_hdf5(self, title: str, pbar: tqdm):
@@ -127,6 +143,7 @@ class MPIGeometryObjects(GeometryObjects[MPIGridView]):
             G: FDTDGrid class describing a grid in a model.
             pbar: Progress bar class instance.
         """
+        assert isinstance(self.grid_view, self.GRID_VIEW_TYPE)
 
         self.grid_view.initialise_materials(self.comm)
 
@@ -142,9 +159,9 @@ class MPIGeometryObjects(GeometryObjects[MPIGridView]):
             self.write_metadata(fdata, title)
 
             dset_slice = (
-                self.grid_view.get_slice(0),
-                self.grid_view.get_slice(1),
-                self.grid_view.get_slice(2),
+                self.grid_view.get_output_slice(0),
+                self.grid_view.get_output_slice(1),
+                self.grid_view.get_output_slice(2),
             )
 
             dset = fdata.create_dataset("/data", self.grid_view.global_size, dtype=data.dtype)
@@ -163,9 +180,9 @@ class MPIGeometryObjects(GeometryObjects[MPIGridView]):
             pbar.update(self.rigidsize)
 
             dset_slice = (
-                self.grid_view.get_slice(0, upper_bound_exclusive=False),
-                self.grid_view.get_slice(1, upper_bound_exclusive=False),
-                self.grid_view.get_slice(2, upper_bound_exclusive=False),
+                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 = fdata.create_dataset(

gprMax/output_controllers/grid_view.py

@@ -54,7 +54,7 @@ class GridView(Generic[GridType]):
 
         self._ID = None
 
-    def get_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice[int, int, int]:
+    def get_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
         """Create a slice object for the specified dimension.
 
         Args:
@@ -272,8 +272,7 @@ class MPIGridView(GridView[MPIGrid]):
         self.global_size = self.size
 
         # Calculate start for the local grid
-        self.global_start = self.start
-        self.start = self.grid.global_to_local_coordinate(self.start)
+        self.global_start = self.grid.local_to_global_coordinate(self.start)
 
         # Bring start into the local grid (and not in the negative halo)
         # local_start must still be aligned with the provided step.
@@ -285,7 +284,7 @@ class MPIGridView(GridView[MPIGrid]):
         )
 
         # Calculate stop for the local grid
-        self.stop = self.grid.global_to_local_coordinate(self.stop)
+        self.global_stop = self.grid.local_to_global_coordinate(self.stop)
 
         self.has_positive_neighbour = self.stop > self.grid.size
 
@@ -309,7 +308,7 @@ class MPIGridView(GridView[MPIGrid]):
         # Update local size
         self.size = self.stop - self.start
 
-    def get_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice[int, int, int]:
+    def get_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
         """Create a slice object for the specified dimension.
 
         Args:
@@ -332,7 +331,19 @@ class MPIGridView(GridView[MPIGrid]):
 
         return slice(self.start[dimension], stop, self.step[dimension])
 
-    def initialise_materials(self, comm: MPI.Comm):
+    def get_output_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
+        if upper_bound_exclusive or self.has_positive_neighbour[dimension]:
+            size = self.size[dimension]
+        else:
+            # Make slice of array one step larger if this rank does not
+            # have a positive neighbour
+            size = self.size[dimension] + 1
+
+        offset = self.offset[dimension]
+
+        return slice(offset, offset + size)
+
+    def initialise_materials(self, comm: MPI.Cartcomm):
         """Create a new ID map for materials in the grid view.
 
         Rather than using the default material IDs (as per the main grid
@@ -349,6 +360,8 @@ class MPIGridView(GridView[MPIGrid]):
 
         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)
@@ -356,25 +369,28 @@ class MPIGridView(GridView[MPIGrid]):
         if materials_by_rank is not None:
             # Filter out duplicate materials and sort by material ID
             all_materials = np.fromiter(chain.from_iterable(materials_by_rank), dtype=Material)
-            unique_materials = np.unique(all_materials)
+            self.materials = np.unique(all_materials)
 
             # The new material IDs corespond to each material's index in
-            # the sorted unique_materials array. For each rank, get the
+            # 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):
-                new_material_ids = np.where(np.isin(unique_materials, materials_by_rank[rank]))[0]
-                comm.Isend([new_material_ids, MPI.INT], rank)
+                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)
 
-            new_material_ids = np.where(np.isin(unique_materials, materials_by_rank[0]))[0]
+            new_material_ids = np.where(np.isin(self.materials, materials_by_rank[0]))[0]
+            new_material_ids = new_material_ids.astype(np.int32)
         else:
-            unique_materials = None
+            self.materials = None
 
             # 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)
 
         # Create map from local material ID to global material ID
-        materials_map = {index: new_id for index, new_id in enumerate(new_material_ids)}
+        materials_map = {
+            local_material_ids[index]: new_id for index, new_id in enumerate(new_material_ids)
+        }
 
         # Create map from material ID to 0 - number of materials
         self.map_materials_func = np.vectorize(lambda id: materials_map[id])

gprMax/user_inputs.py

@@ -180,6 +180,14 @@ class MainGridUserInput(UserInput[GridType]):
 
         return lower_within_grid and upper_within_grid, lower_point, upper_point
 
+    def check_output_object_bounds(
+        self, p1: Tuple[float, float, float], p2: Tuple[float, float, float], cmd_str: str
+    ) -> Tuple[npt.NDArray[np.int32], npt.NDArray[np.int32]]:
+        # We only care if the bounds are in the global grid (an error
+        # will be thrown if that is not the case).
+        _, lower_bound, upper_bound = self._check_2d_points(p1, p2, cmd_str)
+        return lower_bound, upper_bound
+
     def check_box_points(
         self, p1: Tuple[float, float, float], p2: Tuple[float, float, float], cmd_str: str
     ) -> Tuple[bool, npt.NDArray[np.int32], npt.NDArray[np.int32]]:

gprMax/user_objects/cmds_output.py

@@ -1,9 +1,7 @@
 import logging
+from typing import Tuple
 
-from gprMax.geometry_outputs import GeometryObjects as GeometryObjectsUser
-from gprMax.geometry_outputs import MPIGeometryObjects as MPIGeometryObjectsUser
 from gprMax.grid.fdtd_grid import FDTDGrid
-from gprMax.grid.mpi_grid import MPIGrid
 from gprMax.model import Model
 from gprMax.subgrids.grid import SubGridBaseGrid
 from gprMax.user_objects.user_objects import OutputUserObject
@@ -140,41 +138,35 @@ class GeometryObjectsWrite(OutputUserObject):
     def hash(self):
         return "#geometry_objects_write"
 
-    def __init__(self, **kwargs):
-        super().__init__(**kwargs)
+    def __init__(
+        self, p1: Tuple[float, float, float], p2: Tuple[float, float, float], filename: str
+    ):
+        super().__init__(p1=p1, p2=p2, filename=filename)
+        self.lower_bound = p1
+        self.upper_bound = p2
+        self.basefilename = filename
 
     def build(self, model: Model, grid: FDTDGrid):
         if isinstance(grid, SubGridBaseGrid):
-            logger.exception(f"{self.params_str()} do not add geometry objects to subgrids.")
-            raise ValueError
-        try:
-            p1 = self.kwargs["p1"]
-            p2 = self.kwargs["p2"]
-            basefilename = self.kwargs["filename"]
-        except KeyError:
-            logger.exception(f"{self.params_str()} requires exactly seven parameters.")
-            raise
+            raise ValueError(f"{self.params_str()} do not add geometry objects to subgrids.")
 
         uip = self._create_uip(grid)
-        p1, p2 = uip.check_box_points(p1, p2, self.params_str())
-        x0, y0, z0 = p1
-        x1, y1, z1 = p2
 
-        # TODO: Remove these when add parallel build
-        if isinstance(grid, MPIGrid):
-            geometry_object_type = MPIGeometryObjectsUser
-        else:
-            geometry_object_type = GeometryObjectsUser
-
-        g = geometry_object_type(x0, y0, z0, x1, y1, z1, basefilename)
-
-        logger.info(
-            f"Geometry objects in the volume from {p1[0] * grid.dx:g}m, "
-            f"{p1[1] * grid.dy:g}m, {p1[2] * grid.dz:g}m, to "
-            f"{p2[0] * grid.dx:g}m, {p2[1] * grid.dy:g}m, "
-            f"{p2[2] * grid.dz:g}m, will be written to "
-            f"{g.filename_hdf5}, with materials written to "
-            f"{g.filename_materials}"
+        discretised_lower_bound, discretised_upper_bound = uip.check_output_object_bounds(
+            self.lower_bound, self.upper_bound, self.params_str()
         )
 
-        model.geometryobjects.append(g)
+        g = model.add_geometry_object(
+            grid, discretised_lower_bound, discretised_upper_bound, self.basefilename
+        )
+
+        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)
+
+            logger.info(
+                f"Geometry objects in the volume from {p1[0]:g}m,"
+                f" {p1[1]:g}m, {p1[2]:g}m, to {p2[0]:g}m, {p2[1]:g}m,"
+                f" {p2[2]:g}m, will be written to {g.filename_hdf5},"
+                f" with materials written to {g.filename_materials}"
+            )