Update snapshots for parallel build

- Snapshots now use the GridView class
- Snashot type hints have been removed from FDTDGrid and MPIGrid as this
  would cause circular imports

gprMax/__init__.py

@@ -40,13 +40,12 @@ from .user_objects.cmds_multiuse import (
     MaterialRange,
     Rx,
     RxArray,
-    Snapshot,
     SoilPeplinski,
     TransmissionLine,
     VoltageSource,
     Waveform,
 )
-from .user_objects.cmds_output import GeometryObjectsWrite, GeometryView
+from .user_objects.cmds_output import GeometryObjectsWrite, GeometryView, Snapshot
 from .user_objects.cmds_singleuse import (
     Discretisation,
     Domain,

gprMax/grid/fdtd_grid.py

@@ -36,7 +36,6 @@ from gprMax.fractals import FractalVolume
 from gprMax.materials import ListMaterial, Material, PeplinskiSoil, RangeMaterial, process_materials
 from gprMax.pml import CFS, PML, print_pml_info
 from gprMax.receivers import Rx
-from gprMax.snapshots import Snapshot
 from gprMax.sources import HertzianDipole, MagneticDipole, Source, TransmissionLine, VoltageSource
 from gprMax.utilities.utilities import fft_power, get_terminal_width, round_value
 from gprMax.waveforms import Waveform
@@ -111,7 +110,7 @@ class FDTDGrid:
         self.magneticdipoles: List[MagneticDipole] = []
         self.transmissionlines: List[TransmissionLine] = []
         self.rxs: List[Rx] = []
-        self.snapshots: List[Snapshot] = []
+        self.snapshots = []  # List[Snapshot]
 
         self.averagevolumeobjects = True
 

gprMax/grid/mpi_grid.py

@@ -31,7 +31,6 @@ from gprMax.cython.pml_build import pml_sum_er_mr
 from gprMax.grid.fdtd_grid import FDTDGrid
 from gprMax.pml import MPIPML, PML
 from gprMax.receivers import Rx
-from gprMax.snapshots import MPISnapshot, Snapshot
 from gprMax.sources import Source
 
 logger = logging.getLogger(__name__)
@@ -378,7 +377,7 @@ class MPIGrid(FDTDGrid):
         local grid.
         """
         if self.is_coordinator():
-            snapshots_by_rank: List[List[Optional[Snapshot]]] = [[] for _ in range(self.comm.size)]
+            snapshots_by_rank = [[] for _ in range(self.comm.size)]
             for snapshot in self.snapshots:
                 ranks = self.get_ranks_between_coordinates(snapshot.start, snapshot.stop)
                 for rank in range(
@@ -395,9 +394,7 @@ class MPIGrid(FDTDGrid):
         else:
             snapshots_by_rank = None
 
-        snapshots: List[Optional[MPISnapshot]] = self.comm.scatter(
-            snapshots_by_rank, root=self.COORDINATOR_RANK
-        )
+        snapshots = self.comm.scatter(snapshots_by_rank, root=self.COORDINATOR_RANK)
 
         for snapshot in snapshots:
             if snapshot is None:
@@ -507,7 +504,8 @@ class MPIGrid(FDTDGrid):
         Global properties/objects are broadcast to all ranks whereas
         local properties/objects are scattered to the relevant ranks.
         """
-        self.scatter_snapshots()
+        pass
+        # self.scatter_snapshots()
 
         # self._halo_swap_array(self.ID[0])
         # self._halo_swap_array(self.ID[1])

gprMax/hash_cmds_multiuse.py

@@ -31,13 +31,12 @@ from .user_objects.cmds_multiuse import (
     MaterialRange,
     Rx,
     RxArray,
-    Snapshot,
     SoilPeplinski,
     TransmissionLine,
     VoltageSource,
     Waveform,
 )
-from .user_objects.cmds_output import GeometryObjectsWrite, GeometryView
+from .user_objects.cmds_output import GeometryObjectsWrite, GeometryView, Snapshot
 
 logger = logging.getLogger(__name__)
 

gprMax/model.py

@@ -19,7 +19,7 @@
 import datetime
 import logging
 import sys
-from typing import List, Optional, Sequence
+from typing import Dict, List, Optional, Sequence
 
 import humanize
 import numpy as np
@@ -40,12 +40,11 @@ init()
 from tqdm import tqdm
 
 import gprMax.config as config
-
-from .fields_outputs import write_hdf5_outputfile
-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
-from .utilities.utilities import get_terminal_width
+from gprMax.fields_outputs import write_hdf5_outputfile
+from gprMax.grid.fdtd_grid import FDTDGrid
+from gprMax.snapshots import Snapshot, save_snapshots
+from gprMax.utilities.host_info import mem_check_build_all, mem_check_run_all, set_omp_threads
+from gprMax.utilities.utilities import get_terminal_width
 
 logger = logging.getLogger(__name__)
 
@@ -192,8 +191,8 @@ class 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).
+            start: Lower extent of the geometry object (x, y, z).
+            stop: Upper extent of the geometry object (x, y, z).
             basefilename: Output filename of the geometry object.
 
         Returns:
@@ -217,9 +216,9 @@ class 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).
+            start: Lower extent of the geometry view (x, y, z).
+            stop: Upper extent of the geometry view (x, y, z).
+            dl: Discritisation of the geometry view (x, y, z).
             filename: Output filename of the geometry view.
 
         Returns:
@@ -252,8 +251,8 @@ class 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).
+            start: Lower extent of the geometry view (x, y, z).
+            stop: Upper extent of the geometry view (x, y, z).
             filename: Output filename of the geometry view.
 
         Returns:
@@ -272,6 +271,52 @@ class Model:
         self.geometryviews.append(geometry_view)
         return geometry_view
 
+    def add_snapshot(
+        self,
+        grid: FDTDGrid,
+        start: npt.NDArray[np.int32],
+        stop: npt.NDArray[np.int32],
+        dl: npt.NDArray[np.int32],
+        time: int,
+        filename: str,
+        fileext: str,
+        outputs: Dict[str, bool],
+    ) -> Optional[Snapshot]:
+        """Add a snapshot to the provided grid.
+
+        Args:
+            grid: Grid to create a snapshot for.
+            start: Lower extent of the snapshot (x, y, z).
+            stop: Upper extent of the snapshot (x, y, z).
+            dl: Discritisation of the snapshot (x, y, z).
+            time: Iteration number to take the snapshot on
+            filename: Output filename of the snapshot.
+            fileext: File extension of the snapshot.
+            outputs: Fields to use in the snapshot.
+
+        Returns:
+            snapshot: The created snapshot.
+        """
+        snapshot = Snapshot(
+            start[0],
+            start[1],
+            start[2],
+            stop[0],
+            stop[1],
+            stop[2],
+            dl[0],
+            dl[1],
+            dl[2],
+            time,
+            filename,
+            fileext,
+            outputs,
+            grid,
+        )
+        # TODO: Move snapshots into the Model
+        grid.snapshots.append(snapshot)
+        return snapshot
+
     def build(self):
         """Builds the Yee cells for a model."""
 

gprMax/mpi_model.py

@@ -1,5 +1,5 @@
 import logging
-from typing import Optional
+from typing import Dict, Optional
 
 import numpy as np
 import numpy.typing as npt
@@ -12,7 +12,7 @@ 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
+from gprMax.snapshots import MPISnapshot, Snapshot, save_snapshots
 
 logger = logging.getLogger(__name__)
 
@@ -73,8 +73,8 @@ class MPIModel(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)
+            start: Lower extent of the geometry object (x, y, z)
+            stop: Upper extent of the geometry object (x, y, z)
             basefilename: Output filename of the geometry object
 
         Returns:
@@ -107,9 +107,9 @@ class MPIModel(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).
+            start: Lower extent of the geometry view (x, y, z).
+            stop: Upper extent of the geometry view (x, y, z).
+            dl: Discritisation of the geometry view (x, y, z).
             filename: Output filename of the geometry view.
 
         Returns:
@@ -151,8 +151,8 @@ class MPIModel(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).
+            start: Lower extent of the geometry view (x, y, z).
+            stop: Upper extent of the geometry view (x, y, z).
             filename: Output filename of the geometry view.
 
         Returns:
@@ -180,6 +180,60 @@ class MPIModel(Model):
             grid.comm.Split(MPI.UNDEFINED)
             return None
 
+    def add_snapshot(
+        self,
+        grid: MPIGrid,
+        start: npt.NDArray[np.int32],
+        stop: npt.NDArray[np.int32],
+        dl: npt.NDArray[np.int32],
+        time: int,
+        filename: str,
+        fileext: str,
+        outputs: Dict[str, bool],
+    ) -> Optional[MPISnapshot]:
+        """Add a snapshot to the provided grid.
+
+        Args:
+            grid: Grid to create a snapshot for.
+            start: Lower extent of the snapshot (x, y, z).
+            stop: Upper extent of the snapshot (x, y, z).
+            dl: Discritisation of the snapshot (x, y, z).
+            time: Iteration number to take the snapshot on
+            filename: Output filename of the snapshot.
+            fileext: File extension of the snapshot.
+            outputs: Fields to use in the snapshot.
+
+        Returns:
+            snapshot: The new snapshot or None if no snapshot was
+                created.
+        """
+        if grid.local_bounds_overlap_grid(start, stop):
+            snapshot = MPISnapshot(
+                start[0],
+                start[1],
+                start[2],
+                stop[0],
+                stop[1],
+                stop[2],
+                dl[0],
+                dl[1],
+                dl[2],
+                time,
+                filename,
+                fileext,
+                outputs,
+                grid,
+            )
+            # TODO: Move snapshots into the Model
+            grid.snapshots.append(snapshot)
+            return snapshot
+        else:
+            # The MPIGridView created by MPISnapshot will create a new
+            # communicator using MPI_Split. Calling this here prevents
+            # deadlock if not all ranks create the new MPISnapshot.
+            grid.comm.Split(MPI.UNDEFINED)
+            return None
+
     def write_output_data(self):
         """Writes output data, i.e. field data for receivers and snapshots to
         file(s).

gprMax/snapshots.py

@@ -20,16 +20,17 @@ import logging
 import sys
 from enum import IntEnum, unique
 from pathlib import Path
-from typing import Dict, List
+from typing import Dict, Generic, List
 
 import h5py
 import numpy as np
-import numpy.typing as npt
 from evtk.hl import imageToVTK
 from mpi4py import MPI
 from tqdm import tqdm
 
 import gprMax.config as config
+from gprMax.grid.mpi_grid import MPIGrid
+from gprMax.output_controllers.grid_view import GridType, GridView, MPIGridView
 
 from ._version import __version__
 from .cython.snapshots import calculate_snapshot_fields
@@ -52,7 +53,7 @@ def save_snapshots(snapshots: List["Snapshot"]):
     logger.info(f"Snapshot directory: {snapshotdir.resolve()}")
 
     for i, snap in enumerate(snapshots):
-        fn = snapshotdir / Path(snap.filename)
+        fn = snapshotdir / snap.filename
         snap.filename = fn.with_suffix(snap.fileext)
         pbar = tqdm(
             total=snap.nbytes,
@@ -69,7 +70,7 @@ def save_snapshots(snapshots: List["Snapshot"]):
     logger.info("")
 
 
-class Snapshot:
+class Snapshot(Generic[GridType]):
     """Snapshots of the electric and magnetic field values."""
 
     allowableoutputs = {
@@ -95,6 +96,10 @@ class Snapshot:
     # GPU - blocks per grid - set according to largest requested snapshot
     bpg = None
 
+    @property
+    def GRID_VIEW_TYPE(self) -> type[GridView]:
+        return GridView
+
     def __init__(
         self,
         xs: int,
@@ -110,8 +115,7 @@ class Snapshot:
         filename: str,
         fileext: str,
         outputs: Dict[str, bool],
-        grid_dl: npt.NDArray[np.float64],
-        grid_dt: float,
+        grid: GridType,
     ):
         """
         Args:
@@ -124,89 +128,62 @@ class Snapshot:
         """
 
         self.fileext = fileext
-        self.filename = filename
+        self.filename = Path(filename)
         self.time = time
         self.outputs = outputs
-        self.grid_dl = grid_dl
-        self.grid_dt = grid_dt
-
-        self.start = np.array([xs, ys, zs], dtype=np.int32)
-        self.stop = np.array([xf, yf, zf], dtype=np.int32)
-        self.step = np.array([dx, dy, dz], dtype=np.int32)
-        self.size = np.ceil((self.stop - self.start) / self.step).astype(np.int32)
-        self.slice: list[slice] = list(map(slice, self.start, self.stop + self.step, self.step))
+        self.grid_view = self.GRID_VIEW_TYPE(grid, xs, ys, zs, xf, yf, zf, dx, dy, dz)
 
         self.nbytes = 0
 
         # Create arrays to hold the field data for snapshot
         self.snapfields = {}
 
+    @property
+    def grid(self) -> GridType:
+        return self.grid_view.grid
+
     # Properties for backwards compatibility
     @property
     def xs(self) -> int:
-        return self.start[0]
+        return self.grid_view.xs
 
     @property
     def ys(self) -> int:
-        return self.start[1]
+        return self.grid_view.ys
 
     @property
     def zs(self) -> int:
-        return self.start[2]
+        return self.grid_view.zs
 
     @property
     def xf(self) -> int:
-        return self.stop[0]
+        return self.grid_view.xf
 
     @property
     def yf(self) -> int:
-        return self.stop[1]
+        return self.grid_view.yf
 
     @property
     def zf(self) -> int:
-        return self.stop[2]
-
-    @property
-    def dx(self) -> int:
-        return self.step[0]
-
-    @property
-    def dy(self) -> int:
-        return self.step[1]
-
-    @property
-    def dz(self) -> int:
-        return self.step[2]
+        return self.grid_view.zf
 
     @property
     def nx(self) -> int:
-        return self.size[0]
+        return self.grid_view.nx
 
     @property
     def ny(self) -> int:
-        return self.size[1]
+        return self.grid_view.ny
 
     @property
     def nz(self) -> int:
-        return self.size[2]
-
-    @property
-    def sx(self) -> slice:
-        return self.slice[0]
-
-    @property
-    def sy(self) -> slice:
-        return self.slice[1]
-
-    @property
-    def sz(self) -> slice:
-        return self.slice[2]
+        return self.grid_view.nz
 
     def initialise_snapfields(self):
         for k, v in self.outputs.items():
             if v:
                 self.snapfields[k] = np.zeros(
-                    (self.nx, self.ny, self.nz),
+                    self.grid_view.size,
                     dtype=config.sim_config.dtypes["float_or_double"],
                 )
                 self.nbytes += self.snapfields[k].nbytes
@@ -217,7 +194,7 @@ class Snapshot:
                     (1, 1, 1), dtype=config.sim_config.dtypes["float_or_double"]
                 )
 
-    def store(self, G):
+    def store(self):
         """Store (in memory) electric and magnetic field values for snapshot.
 
         Args:
@@ -225,12 +202,12 @@ class Snapshot:
         """
 
         # Memory views of field arrays to dimensions required for the snapshot
-        Exslice = np.ascontiguousarray(G.Ex[self.sx, self.sy, self.sz])
-        Eyslice = np.ascontiguousarray(G.Ey[self.sx, self.sy, self.sz])
-        Ezslice = np.ascontiguousarray(G.Ez[self.sx, self.sy, self.sz])
-        Hxslice = np.ascontiguousarray(G.Hx[self.sx, self.sy, self.sz])
-        Hyslice = np.ascontiguousarray(G.Hy[self.sx, self.sy, self.sz])
-        Hzslice = np.ascontiguousarray(G.Hz[self.sx, self.sy, self.sz])
+        Exslice = self.grid_view.get_Ex()
+        Eyslice = self.grid_view.get_Ey()
+        Ezslice = self.grid_view.get_Ez()
+        Hxslice = self.grid_view.get_Hx()
+        Hyslice = self.grid_view.get_Hy()
+        Hzslice = self.grid_view.get_Hz()
 
         # Calculate field values at points (comes from averaging field
         # components in cells)
@@ -286,10 +263,13 @@ class Snapshot:
             if self.outputs.get(k)
         }
 
+        origin = self.grid_view.start * self.grid.dl
+        spacing = self.grid_view.step * self.grid.dl
+
         imageToVTK(
             str(self.filename.with_suffix("")),
-            origin=tuple(self.start * self.step * self.grid_dl),
-            spacing=tuple(self.step * self.grid_dl),
+            origin=tuple(origin),
+            spacing=tuple(spacing),
             cellData=celldata,
         )
 
@@ -312,9 +292,9 @@ class Snapshot:
         f.attrs["gprMax"] = __version__
         # TODO: Output model name (title) and grid name? in snapshot output
         # f.attrs["Title"] = G.title
-        f.attrs["nx_ny_nz"] = (self.nx, self.ny, self.nz)
-        f.attrs["dx_dy_dz"] = self.step * self.grid_dl
-        f.attrs["time"] = self.time * self.grid_dt
+        f.attrs["nx_ny_nz"] = tuple(self.grid_view.size)
+        f.attrs["dx_dy_dz"] = self.grid_view.step * self.grid.dl
+        f.attrs["time"] = self.time * self.grid.dt
 
         for key in ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"]:
             if self.outputs[key]:
@@ -337,12 +317,16 @@ class Dir(IntEnum):
     POS = 1
 
 
-class MPISnapshot(Snapshot):
+class MPISnapshot(Snapshot[MPIGrid]):
     H_TAG = 0
     EX_TAG = 1
     EY_TAG = 2
     EZ_TAG = 3
 
+    @property
+    def GRID_VIEW_TYPE(self) -> type[MPIGridView]:
+        return MPIGridView
+
     def __init__(
         self,
         xs: int,
@@ -358,27 +342,23 @@ class MPISnapshot(Snapshot):
         filename: str,
         fileext: str,
         outputs: Dict[str, bool],
-        grid_dl: npt.NDArray[np.float64],
-        grid_dt: float,
+        grid: MPIGrid,
     ):
-        super().__init__(
-            xs, ys, zs, xf, yf, zf, dx, dy, dz, time, filename, fileext, outputs, grid_dl, grid_dt
-        )
+        super().__init__(xs, ys, zs, xf, yf, zf, dx, dy, dz, time, filename, fileext, outputs, grid)
 
-        self.offset = np.zeros(3, dtype=np.int32)
-        self.global_size = self.size.copy()
+        assert isinstance(self.grid_view, self.GRID_VIEW_TYPE)
+        self.comm = self.grid_view.comm
 
-        self.comm: MPI.Cartcomm = None  # type: ignore
-
-    def initialise_snapfields(self):
-        # Start and stop may have changed since initialisation
-        self.size = np.ceil((self.stop - self.start) / self.step).astype(np.int32)
-        return super().initialise_snapfields()
+        # Get neighbours
+        self.neighbours = np.full((3, 2), -1, dtype=int)
+        self.neighbours[Dim.X] = self.comm.Shift(direction=Dim.X, disp=1)
+        self.neighbours[Dim.Y] = self.comm.Shift(direction=Dim.Y, disp=1)
+        self.neighbours[Dim.Z] = self.comm.Shift(direction=Dim.Z, disp=1)
 
     def has_neighbour(self, dimension: Dim, direction: Dir) -> bool:
         return self.neighbours[dimension][direction] != -1
 
-    def store(self, G):
+    def store(self):
         """Store (in memory) electric and magnetic field values for snapshot.
 
         Args:
@@ -387,31 +367,17 @@ class MPISnapshot(Snapshot):
 
         logger.debug(f"Saving snapshot for iteration: {self.time}")
 
-        # Get neighbours
-        self.neighbours = np.full((3, 2), -1, dtype=int)
-        self.neighbours[Dim.X] = self.comm.Shift(direction=Dim.X, disp=1)
-        self.neighbours[Dim.Y] = self.comm.Shift(direction=Dim.Y, disp=1)
-        self.neighbours[Dim.Z] = self.comm.Shift(direction=Dim.Z, disp=1)
-
-        # If we do not have a positive neighbour, add an extra step to
-        # make the upper bound inclusive. Otherwise the additional step
-        # will be provided by the received halo.
-        slice_stop = np.where(
-            self.neighbours[:, Dir.POS] == -1,
-            self.stop + self.step,
-            self.stop,
-        )
-        self.slice = list(map(slice, self.start, slice_stop, self.step))
-
         # Memory views of field arrays to dimensions required for the snapshot
-        Exslice = np.ascontiguousarray(G.Ex[self.sx, self.sy, self.sz])
-        Eyslice = np.ascontiguousarray(G.Ey[self.sx, self.sy, self.sz])
-        Ezslice = np.ascontiguousarray(G.Ez[self.sx, self.sy, self.sz])
-        Hxslice = np.ascontiguousarray(G.Hx[self.sx, self.sy, self.sz])
-        Hyslice = np.ascontiguousarray(G.Hy[self.sx, self.sy, self.sz])
-        Hzslice = np.ascontiguousarray(G.Hz[self.sx, self.sy, self.sz])
+        Exslice = self.grid_view.get_Ex()
+        Eyslice = self.grid_view.get_Ey()
+        Ezslice = self.grid_view.get_Ez()
+        Hxslice = self.grid_view.get_Hx()
+        Hyslice = self.grid_view.get_Hy()
+        Hzslice = self.grid_view.get_Hz()
 
         """
+        Halos required by each field to average field components:
+
         Exslice - y + z halo
         Eyslice - x + z halo
         Ezslice - x + y halo
@@ -579,25 +545,23 @@ class MPISnapshot(Snapshot):
         Args:
             pbar: Progress bar class instance.
         """
+        assert isinstance(self.grid_view, self.GRID_VIEW_TYPE)
 
         f = h5py.File(self.filename, "w", driver="mpio", comm=self.comm)
 
         f.attrs["gprMax"] = __version__
         # TODO: Output model name (title) and grid name? in snapshot output
         # f.attrs["Title"] = G.title
-        f.attrs["nx_ny_nz"] = self.global_size
-        f.attrs["dx_dy_dz"] = self.step * self.grid_dl
-        f.attrs["time"] = self.time * self.grid_dt
+        f.attrs["nx_ny_nz"] = self.grid_view.global_size
+        f.attrs["dx_dy_dz"] = self.grid_view.step * self.grid.dl
+        f.attrs["time"] = self.time * self.grid.dt
+
+        dset_slice = self.grid_view.get_3d_output_slice()
 
         for key in ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"]:
             if self.outputs[key]:
-                dset = f.create_dataset(key, self.global_size)
-                # TODO: Is there a better way to do this slice?
-                start = self.offset
-                stop = start + self.size
-                dset[start[0] : stop[0], start[1] : stop[1], start[2] : stop[2]] = self.snapfields[
-                    key
-                ]
+                dset = f.create_dataset(key, self.grid_view.global_size)
+                dset[dset_slice] = self.snapfields[key]
                 pbar.update(n=self.snapfields[key].nbytes)
 
         f.close()

gprMax/updates/cpu_updates.py

@@ -50,7 +50,7 @@ class CPUUpdates(Updates[GridType]):
         """
         for snap in self.grid.snapshots:
             if snap.time == iteration + 1:
-                snap.store(self.grid)
+                snap.store()
 
     def update_magnetic(self):
         """Updates magnetic field components."""

gprMax/user_inputs.py

@@ -50,15 +50,10 @@ class UserInput(Generic[GridType]):
     def __init__(self, grid: GridType):
         self.grid = grid
 
-    def point_within_bounds(
-        self, p: npt.NDArray[np.int32], cmd_str: str, name: str = "", ignore_error=False
-    ) -> bool:
+    def point_within_bounds(self, p: npt.NDArray[np.int32], cmd_str: str, name: str = "") -> bool:
         try:
             return self.grid.within_bounds(p)
         except ValueError as err:
-            if ignore_error:
-                return False
-
             v = ["x", "y", "z"]
             # Discretisation
             dl = getattr(self.grid, f"d{err.args[0]}")
@@ -71,9 +66,6 @@ class UserInput(Generic[GridType]):
             logger.exception(s)
             raise
 
-    def grid_upper_bound(self) -> list[int]:
-        return [self.grid.nx, self.grid.ny, self.grid.nz]
-
     def discretise_static_point(self, point: Tuple[float, float, float]) -> npt.NDArray[np.int32]:
         """Get the nearest grid index to a continuous static point.
 

gprMax/user_objects/cmds_multiuse.py

@@ -1054,217 +1054,6 @@ class RxArray(GridUserObject):
                     receiver.build(grid)
 
 
-class Snapshot(GridUserObject):
-    """Obtains information about the electromagnetic fields within a volume
-        of the model at a given time instant.
-
-    Attributes:
-        p1: tuple required to specify lower left (x,y,z) coordinates of volume
-                of snapshot in metres.
-        p2: tuple required to specify upper right (x,y,z) coordinates of volume
-                of snapshot in metres.
-        dl: tuple require to specify spatial discretisation of the snapshot
-                in metres.
-        filename: string required for name of the file to store snapshot.
-        time/iterations: either a float for time or an int for iterations
-                            must be specified for point in time at which the
-                            snapshot will be taken.
-        fileext: optional string to indicate type for snapshot file, either
-                            '.vti' (default) or '.h5'
-        outputs: optional list of outputs for receiver. It can be any
-                    selection from Ex, Ey, Ez, Hx, Hy, or Hz.
-    """
-
-    # TODO: Make this an output user object
-    @property
-    def order(self):
-        return 9
-
-    @property
-    def hash(self):
-        return "#snapshot"
-
-    def __init__(self, **kwargs):
-        super().__init__(**kwargs)
-
-    def _calculate_upper_bound(
-        self, start: npt.NDArray[np.int32], step: npt.NDArray[np.int32], size: npt.NDArray[np.int32]
-    ) -> npt.NDArray[np.int32]:
-        # upper_bound = p2 + dl - ((snapshot_size - 1) % dl) - 1
-        return start + step * np.ceil(size / step)
-
-    def build(self, grid: FDTDGrid):
-        if isinstance(grid, SubGridBaseGrid):
-            logger.exception(f"{self.params_str()} do not add snapshots to subgrids.")
-            raise ValueError
-        try:
-            p1 = self.kwargs["p1"]
-            p2 = self.kwargs["p2"]
-            dl = self.kwargs["dl"]
-            filename = self.kwargs["filename"]
-        except KeyError:
-            logger.exception(f"{self.params_str()} requires exactly 11 parameters.")
-            raise
-
-        uip = self._create_uip(grid)
-        dl = uip.discretise_static_point(dl)
-
-        try:
-            _, 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
-
-        snapshot_size = p2 - p1
-
-        # If p2 does not line up with the set discretisation, the actual
-        # maximum element accessed in the grid will be this upper bound.
-        upper_bound = self._calculate_upper_bound(p1, dl, snapshot_size)
-
-        # Each coordinate may need a different method to correct p2.
-        # Therefore, this check needs to be repeated after each
-        # correction has been applied.
-        while any(p2 < upper_bound):
-            # Ideally extend p2 up to the correct upper bound. This will
-            # not change the snapshot output.
-            if uip.point_within_bounds(upper_bound, "", "", ignore_error=True):
-                p2 = upper_bound
-                p2_continuous = p2 * grid.dl
-                logger.warning(
-                    f"{self.params_str()} upper bound not aligned with discretisation. Updating 'p2'"
-                    f" to {p2_continuous}"
-                )
-            # If the snapshot size cannot be increased, the
-            # discretisation may need reducing. E.g. for snapshots of 2D
-            # models.
-            elif any(dl > snapshot_size):
-                dl = np.where(dl > snapshot_size, snapshot_size, dl)
-                upper_bound = self._calculate_upper_bound(p1, dl, snapshot_size)
-                dl_continuous = dl * grid.dl
-                logger.warning(
-                    f"{self.params_str()} current bounds and discretisation would go outside"
-                    f" domain. As discretisation is larger than the snapshot size in at least one"
-                    f" dimension, limiting 'dl' to {dl_continuous}"
-                )
-            # Otherwise, limit p2 to the discretisation step below the
-            # current snapshot size. This will reduce the size of the
-            # snapshot by 1 in the effected dimension(s), but avoid out
-            # of memory access.
-            else:
-                p2 = np.where(uip.grid_upper_bound() < upper_bound, p2 - (snapshot_size % dl), p2)
-                snapshot_size = p2 - p1
-                upper_bound = self._calculate_upper_bound(p1, dl, snapshot_size)
-                p2_continuous = p2 * grid.dl
-                logger.warning(
-                    f"{self.params_str()} current bounds and discretisation would go outside"
-                    f" domain. Limiting 'p2' to {p2_continuous}"
-                )
-
-        if any(dl < 0):
-            logger.exception(f"{self.params_str()} the step size should not be less than zero.")
-            raise ValueError
-        if any(dl < 1):
-            logger.exception(
-                f"{self.params_str()} the step size should not be less than the spatial discretisation."
-            )
-            raise ValueError
-
-        # If number of iterations given
-        try:
-            iterations = self.kwargs["iterations"]
-        # If real floating point value given
-        except KeyError:
-            try:
-                time = self.kwargs["time"]
-            except KeyError:
-                logger.exception(f"{self.params_str()} requires exactly 5 parameters.")
-                raise
-            if time > 0:
-                iterations = round_value((time / grid.dt)) + 1
-            else:
-                logger.exception(f"{self.params_str()} time value must be greater than zero.")
-                raise ValueError
-
-        if iterations <= 0 or iterations > grid.iterations:
-            logger.exception(f"{self.params_str()} time value is not valid.")
-            raise ValueError
-
-        try:
-            fileext = self.kwargs["fileext"]
-            if fileext not in SnapshotUser.fileexts:
-                logger.exception(
-                    f"'{fileext}' is not a valid format for a "
-                    "snapshot file. Valid options are: "
-                    f"{' '.join(SnapshotUser.fileexts)}."
-                )
-                raise ValueError
-        except KeyError:
-            fileext = SnapshotUser.fileexts[0]
-
-        if isinstance(grid, MPIGrid) and fileext != ".h5":
-            logger.exception(
-                f"{self.params_str()} currently only '.h5' snapshots are compatible with MPI."
-            )
-            raise ValueError
-
-        try:
-            tmp = self.kwargs["outputs"]
-            outputs = dict.fromkeys(SnapshotUser.allowableoutputs, False)
-            # Check and set output names
-            for output in tmp:
-                if output not in SnapshotUser.allowableoutputs.keys():
-                    logger.exception(
-                        f"{self.params_str()} contains an output "
-                        f"type that is not allowable. Allowable "
-                        f"outputs in current context are "
-                        f"{', '.join(SnapshotUser.allowableoutputs.keys())}."
-                    )
-                    raise ValueError
-                else:
-                    outputs[output] = True
-        except KeyError:
-            # If outputs are not specified, use default
-            outputs = dict.fromkeys(SnapshotUser.allowableoutputs, True)
-
-        if isinstance(grid, MPIGrid):
-            snapshot_type = MPISnapshotUser
-        else:
-            snapshot_type = SnapshotUser
-
-        xs, ys, zs = p1
-        xf, yf, zf = p2
-        dx, dy, dz = dl
-
-        s = snapshot_type(
-            xs,
-            ys,
-            zs,
-            xf,
-            yf,
-            zf,
-            dx,
-            dy,
-            dz,
-            iterations,
-            filename,
-            fileext=fileext,
-            outputs=outputs,
-            grid_dl=grid.dl,
-            grid_dt=grid.dt,
-        )
-
-        logger.info(
-            f"Snapshot from {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to "
-            f"{xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m, discretisation "
-            f"{dx * grid.dx:g}m, {dy * grid.dy:g}m, {dz * grid.dz:g}m, "
-            f"at {s.time * grid.dt:g} secs with field outputs "
-            f"{', '.join([k for k, v in outputs.items() if v])} and "
-            f"filename {s.filename}{s.fileext} will be created."
-        )
-
-        grid.snapshots.append(s)
-
-
 class Material(GridUserObject):
     """Specifies a material in the model described by a set of constitutive
         parameters.

gprMax/user_objects/cmds_output.py

@@ -1,14 +1,249 @@
 import logging
-from typing import Tuple
+from typing import List, Optional, Tuple
+
+import numpy as np
+import numpy.typing as npt
 
 from gprMax.grid.fdtd_grid import FDTDGrid
+from gprMax.grid.mpi_grid import MPIGrid
 from gprMax.model import Model
+from gprMax.snapshots import Snapshot as SnapshotUser
 from gprMax.subgrids.grid import SubGridBaseGrid
 from gprMax.user_objects.user_objects import OutputUserObject
+from gprMax.utilities.utilities import round_int
 
 logger = logging.getLogger(__name__)
 
 
+class Snapshot(OutputUserObject):
+    """Obtains information about the electromagnetic fields within a volume
+        of the model at a given time instant.
+
+    Attributes:
+        p1: tuple required to specify lower left (x,y,z) coordinates of volume
+                of snapshot in metres.
+        p2: tuple required to specify upper right (x,y,z) coordinates of volume
+                of snapshot in metres.
+        dl: tuple require to specify spatial discretisation of the snapshot
+                in metres.
+        filename: string required for name of the file to store snapshot.
+        time/iterations: either a float for time or an int for iterations
+                            must be specified for point in time at which the
+                            snapshot will be taken.
+        fileext: optional string to indicate type for snapshot file, either
+                            '.vti' (default) or '.h5'
+        outputs: optional list of outputs for receiver. It can be any
+                    selection from Ex, Ey, Ez, Hx, Hy, or Hz.
+    """
+
+    @property
+    def order(self):
+        return 9
+
+    @property
+    def hash(self):
+        return "#snapshot"
+
+    def __init__(
+        self,
+        p1: Tuple[float, float, float],
+        p2: Tuple[float, float, float],
+        dl: Tuple[float, float, float],
+        filename: str,
+        fileext: Optional[str] = None,
+        iterations: Optional[int] = None,
+        time: Optional[float] = None,
+        outputs: Optional[List[str]] = None,
+    ):
+        super().__init__(
+            p1=p1,
+            p2=p2,
+            dl=dl,
+            filename=filename,
+            fileext=fileext,
+            iterations=iterations,
+            time=time,
+            outputs=outputs,
+        )
+        self.lower_bound = p1
+        self.upper_bound = p2
+        self.dl = dl
+        self.filename = filename
+        self.file_extension = fileext
+        self.iterations = iterations
+        self.time = time
+        self.outputs = outputs
+
+    def _calculate_upper_bound(
+        self, start: npt.NDArray[np.int32], step: npt.NDArray[np.int32], size: npt.NDArray[np.int32]
+    ) -> npt.NDArray[np.int32]:
+        return start + step * np.ceil(size / step)
+
+    def build(self, model: Model, grid: FDTDGrid):
+        if isinstance(grid, SubGridBaseGrid):
+            raise ValueError(f"{self.params_str()} do not add snapshots to subgrids.")
+
+        uip = self._create_uip(grid)
+        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)
+
+        snapshot_size = discretised_upper_bound - discretised_lower_bound
+
+        # If p2 does not line up with the set discretisation, the actual
+        # maximum element accessed in the grid will be this upper bound.
+        upper_bound = self._calculate_upper_bound(
+            discretised_lower_bound, discretised_dl, snapshot_size
+        )
+
+        # Each coordinate may need a different method to correct p2.
+        # Therefore, this check needs to be repeated after each
+        # correction has been applied.
+        while any(discretised_upper_bound < upper_bound):
+            try:
+                uip.point_within_bounds(
+                    upper_bound, f"[{upper_bound[0]}, {upper_bound[1]}, {upper_bound[2]}]"
+                )
+                upper_bound_within_grid = True
+            except ValueError:
+                upper_bound_within_grid = False
+
+            # Ideally extend p2 up to the correct upper bound. This will
+            # not change the snapshot output.
+            if upper_bound_within_grid:
+                discretised_upper_bound = upper_bound
+                upper_bound_continuous = discretised_upper_bound * grid.dl
+                logger.warning(
+                    f"{self.params_str()} upper bound not aligned with discretisation. Updating 'p2'"
+                    f" to {upper_bound_continuous}"
+                )
+            # If the snapshot size cannot be increased, the
+            # discretisation may need reducing. E.g. for snapshots of 2D
+            # models.
+            elif any(discretised_dl > snapshot_size):
+                discretised_dl = np.where(
+                    discretised_dl > snapshot_size, snapshot_size, discretised_dl
+                )
+                upper_bound = self._calculate_upper_bound(
+                    discretised_lower_bound, discretised_dl, snapshot_size
+                )
+                dl_continuous = discretised_dl * grid.dl
+                logger.warning(
+                    f"{self.params_str()} current bounds and discretisation would go outside"
+                    f" domain. As discretisation is larger than the snapshot size in at least one"
+                    f" dimension, limiting 'dl' to {dl_continuous}"
+                )
+            # Otherwise, limit p2 to the discretisation step below the
+            # current snapshot size. This will reduce the size of the
+            # snapshot by 1 in the effected dimension(s), but avoid out
+            # of memory access.
+            else:
+                discretised_upper_bound = np.where(
+                    discretised_upper_bound < upper_bound,
+                    upper_bound - discretised_dl,
+                    discretised_upper_bound,
+                )
+                snapshot_size = discretised_upper_bound - discretised_lower_bound
+                upper_bound = self._calculate_upper_bound(
+                    discretised_lower_bound, discretised_dl, snapshot_size
+                )
+                upper_bound_continuous = discretised_upper_bound * grid.dl
+                logger.warning(
+                    f"{self.params_str()} current bounds and discretisation would go outside"
+                    f" domain. Limiting 'p2' to {upper_bound_continuous}"
+                )
+
+                # Raise error to prevent an infinite loop. This is here
+                # as a precaution, it shouldn't be needed.
+                if any(discretised_upper_bound < upper_bound):
+                    raise ValueError(f"{self.params_str()} invalid snapshot.")
+
+        if any(discretised_dl < 0):
+            raise ValueError(f"{self.params_str()} the step size should not be less than zero.")
+        if any(discretised_dl < 1):
+            raise ValueError(
+                f"{self.params_str()} the step size should not be less than the spatial discretisation."
+            )
+
+        if self.iterations is not None and self.time is not None:
+            logger.warning(
+                f"{self.params_str()} Time and iterations were both specified, using 'iterations'"
+            )
+
+        # If number of iterations given
+        if self.iterations is not None:
+            if self.iterations <= 0 or self.iterations > grid.iterations:
+                raise ValueError(f"{self.params_str()} time value is not valid.")
+
+        # If time value given
+        elif self.time is not None:
+            if self.time > 0:
+                self.iterations = round_int((self.time / grid.dt)) + 1
+            else:
+                raise ValueError(f"{self.params_str()} time value must be greater than zero.")
+
+        # No iteration or time value given
+        else:
+            raise ValueError(f"{self} specify a time or number of iterations")
+
+        if self.file_extension is None:
+            self.file_extension = SnapshotUser.fileexts[0]
+        elif self.file_extension not in SnapshotUser.fileexts:
+            raise ValueError(
+                f"'{self.file_extension}' is not a valid format for a snapshot file."
+                f" Valid options are: {' '.join(SnapshotUser.fileexts)}."
+            )
+
+        # TODO: Allow VTKHDF files when they are implemented
+        if isinstance(grid, MPIGrid) and self.file_extension != ".h5":
+            raise ValueError(
+                f"{self.params_str()} currently only '.h5' snapshots are compatible with MPI."
+            )
+
+        if self.outputs is None:
+            outputs = dict.fromkeys(SnapshotUser.allowableoutputs, True)
+        else:
+            outputs = dict.fromkeys(SnapshotUser.allowableoutputs, False)
+            # Check and set output names
+            for output in self.outputs:
+                if output not in SnapshotUser.allowableoutputs.keys():
+                    raise ValueError(
+                        f"{self.params_str()} contains an output type that is not"
+                        " allowable. Allowable outputs in current context are "
+                        f"{', '.join(SnapshotUser.allowableoutputs.keys())}."
+                    )
+                else:
+                    outputs[output] = True
+
+        snapshot = model.add_snapshot(
+            grid,
+            discretised_lower_bound,
+            discretised_upper_bound,
+            discretised_dl,
+            self.iterations,
+            self.filename,
+            self.file_extension,
+            outputs,
+        )
+
+        if snapshot is not None:
+            p1 = uip.round_to_grid_static_point(self.lower_bound)
+            p2 = uip.round_to_grid_static_point(self.upper_bound)
+            dl = uip.round_to_grid_static_point(self.dl)
+
+            logger.info(
+                f"{self.grid_name(grid)}Snapshot from"
+                f" {p1[0]:g}m, {p1[1]:g}m, {p1[2]:g}m, to"
+                f" {p2[0]:g}m, {p2[1]:g}m, {p2[2]:g}m, discretisation"
+                f" {dl[0]:g}m, {dl[1]:g}m, {dl[2]:g}m, at"
+                f" {snapshot.time * grid.dt:g} secs with field outputs"
+                f" {', '.join([k for k, v in outputs.items() if v])} "
+                f" and filename {snapshot.filename}{snapshot.fileext}"
+                " will be created."
+            )
+
+
 class GeometryView(OutputUserObject):
     """Outputs to file(s) information about the geometry (mesh) of model.