Add MPIFractalSurface class

gprMax/fractals/fractal_surface.py

@@ -1,10 +1,18 @@
+import logging
+from typing import List, Optional, Tuple
+
 import numpy as np
 import numpy.typing as npt
+from mpi4py import MPI
+from mpi4py_fft import PFFT, DistArray, newDistArray
+from mpi4py_fft.pencil import Subcomm
 from scipy import fftpack
 
 from gprMax import config
 from gprMax.cython.fractals_generate import generate_fractal2D
+from gprMax.utilities.mpi import Dim, Dir, get_relative_neighbour
 
+logger = logging.getLogger(__name__)
 np.seterr(divide="raise")
 
 
@@ -169,3 +177,300 @@ class FractalSurface:
             + self.fractalrange[0]
             - ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * fractalmin
         )
+
+
+class MPIFractalSurface(FractalSurface):
+    def __init__(
+        self,
+        xs: int,
+        xf: int,
+        ys: int,
+        yf: int,
+        zs: int,
+        zf: int,
+        dimension: float,
+        seed: Optional[int],
+        comm: MPI.Cartcomm,
+        ux: int,
+        uy: int,
+        uz: int,
+    ):
+        super().__init__(xs, xf, ys, yf, zs, zf, dimension, seed)
+        self.comm = comm
+        self.upper_bound = np.array([ux, uy, uz])
+
+    def generate_fractal_surface(self):
+        """Generate a 2D array with a fractal distribution."""
+
+        if self.xs == self.xf:
+            color = self.xs
+            static_dimension = Dim.X
+            dims = [Dim.Y, Dim.Z]
+        elif self.ys == self.yf:
+            color = self.ys
+            static_dimension = Dim.Y
+            dims = [Dim.X, Dim.Z]
+        elif self.zs == self.zf:
+            color = self.zs
+            static_dimension = Dim.Z
+            dims = [Dim.X, Dim.Y]
+
+        # Exit early if this rank does not contain the Fractal Surface
+        if (
+            any(self.stop[dims] < 0)
+            or any(self.start[dims] > self.upper_bound[dims])
+            or self.fractalrange[1] < 0
+            or self.fractalrange[0] > self.upper_bound[static_dimension]
+        ):
+            self.comm.Split(MPI.UNDEFINED)
+            # Update start and stop to local bounds
+            self.start = np.maximum(self.start, 0)
+            self.start = np.minimum(self.start, self.upper_bound)
+            self.stop = np.maximum(self.stop, 0)
+            self.stop = np.minimum(self.stop, self.upper_bound)
+            return
+        else:
+            # Create new cartsesian communicator for the Fractal Surface
+            comm = self.comm.Split(color=color)
+            assert isinstance(comm, MPI.Intracomm)
+            min_coord = np.array(self.comm.coords, dtype=np.int32)[dims]
+            max_coord = min_coord + 1
+            comm.Allreduce(MPI.IN_PLACE, min_coord, MPI.MIN)
+            comm.Allreduce(MPI.IN_PLACE, max_coord, MPI.MAX)
+            cart_dims = (max_coord - min_coord).tolist()
+            self.comm = comm.Create_cart(cart_dims)
+
+        # Check domain decomosition is valid for the Fractal Volume
+        if all([dim > 1 for dim in self.comm.dims]):
+            raise ValueError(
+                "Fractal surface must be positioned such that its MPI decomposition is 1 in at least"
+                f" 1 dimension. Current decompostion is: {self.comm.dims}"
+            )
+
+        surfacedims = self.get_surface_dims()
+
+        # Positional vector at centre of array, scaled by weighting
+        v1 = self.weighting * surfacedims / 2
+
+        subcomm = Subcomm(self.comm)
+
+        A = DistArray(self.size[dims], subcomm, dtype=self.dtype)
+
+        fft = PFFT(
+            None,
+            axes=tuple(np.argsort(self.comm.dims)[::-1]),
+            darray=A,
+            collapse=False,
+            backend="fftw",
+        )
+
+        # Decomposition of A may be different to the MPIGrid
+        A_shape = np.array(A.shape)
+        A_substart = np.array(A.substart)
+
+        # 3D array of random numbers to be convolved with the fractal function
+        rng = np.random.default_rng(seed=self.seed)
+
+        for index in np.ndindex(*A.global_shape):
+            index = np.array(index)
+            if any(index < A_substart) or any(index >= A_substart + A_shape):
+                rng.standard_normal()
+            else:
+                index -= A_substart
+                A[index[0], index[1]] = rng.standard_normal()
+
+        A_hat = newDistArray(fft)
+        assert isinstance(A_hat, DistArray)
+
+        # 2D FFT
+        fft.forward(A, A_hat, normalize=False)
+
+        # Generate fractal
+        generate_fractal2D(
+            A_hat.shape[0],
+            A_hat.shape[1],
+            A_hat.substart[0],
+            A_hat.substart[1],
+            A_hat.global_shape[0],
+            A_hat.global_shape[1],
+            config.get_model_config().ompthreads,
+            self.dimension,
+            self.weighting,
+            v1,
+            A_hat,
+            A_hat,
+        )
+
+        # Set DC component of FFT to zero
+        if all(A_substart == 0):
+            A_hat[0, 0] = 0
+
+        # Inverse 2D FFT transform
+        fft.backward(A_hat, A, normalize=True)
+
+        # Take the real part (numerical errors can give rise to an imaginary part)
+        #  of the IFFT, and convert type to floattype. N.B calculation of fractals
+        # must always be carried out at double precision, i.e. float64, complex128
+        A = np.real(A).astype(config.sim_config.dtypes["float_or_double"], copy=False)
+
+        # Allreduce to get min and max values in the fractal surface
+        min_value = np.array(np.amin(A), dtype=config.sim_config.dtypes["float_or_double"])
+        max_value = np.array(np.amax(A), dtype=config.sim_config.dtypes["float_or_double"])
+        self.comm.Allreduce(MPI.IN_PLACE, min_value, MPI.MIN)
+        self.comm.Allreduce(MPI.IN_PLACE, max_value, MPI.MAX)
+
+        # Scale the fractal volume according to requested range
+        fractalrange = max_value - min_value
+        A = (
+            A * ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange)
+            + self.fractalrange[0]
+            - ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * min_value
+        )
+
+        # Distribute A (DistArray) to match the MPIGrid decomposition
+        local_shape = (np.minimum(self.stop, self.upper_bound) - np.maximum(self.start, 0))[dims]
+        self.fractalsurface = np.zeros(
+            local_shape,
+            dtype=config.sim_config.dtypes["float_or_double"],
+        )
+
+        # Negative means send to negative neighbour
+        # Positive means receive from negative neighbour
+        negative_offset = np.where(self.start[dims] >= 0, 0, self.start[dims] + A_substart)
+
+        # Negative means send to positive neighbour
+        # Positive means receive from positive neighbour
+        positive_offset = np.minimum(self.stop, self.upper_bound)[dims] - (
+            self.start[dims] + A_substart + A_shape
+        )
+
+        dirs = np.full(2, Dir.NONE)
+
+        starts, subshape = self.calculate_starts_and_subshape(
+            A_shape, -negative_offset, -positive_offset, dirs, sending=True
+        )
+        ends = starts + subshape
+        A_local = A[starts[0] : ends[0], starts[1] : ends[1]]
+
+        starts, subshape = self.calculate_starts_and_subshape(
+            local_shape, negative_offset, positive_offset, dirs
+        )
+        ends = starts + subshape
+        self.fractalsurface[starts[0] : ends[0], starts[1] : ends[1]] = A_local
+
+        requests: List[MPI.Request] = []
+
+        # Need to check neighbours in each direction (2D plane)
+        sections = [
+            (Dir.NEG, Dir.NONE),
+            (Dir.POS, Dir.NONE),
+            (Dir.NONE, Dir.NEG),
+            (Dir.NONE, Dir.POS),
+            (Dir.NEG, Dir.NEG),
+            (Dir.NEG, Dir.POS),
+            (Dir.POS, Dir.NEG),
+            (Dir.POS, Dir.POS),
+        ]
+
+        for section in sections:
+            dirs[0] = section[0]
+            dirs[1] = section[1]
+            rank = get_relative_neighbour(self.comm, dirs)
+
+            # Skip if no neighbour
+            if rank == -1:
+                continue
+
+            # Check if any data to send
+            if all(
+                np.logical_or(
+                    dirs == Dir.NONE,
+                    np.where(dirs == Dir.NEG, negative_offset <= 0, positive_offset <= 0),
+                )
+            ):
+                mpi_type = self.create_mpi_type(
+                    A_shape, -negative_offset, -positive_offset, dirs, sending=True
+                )
+
+                logger.debug(
+                    f"Sending fractal surface to rank {rank}, MPI type={mpi_type.decode()}"
+                )
+                self.comm.Isend([A, mpi_type], rank)
+
+            # Check if any data to receive
+            if all(
+                np.logical_or(
+                    dirs == Dir.NONE,
+                    np.where(dirs == Dir.NEG, negative_offset > 0, positive_offset > 0),
+                )
+            ):
+                mpi_type = self.create_mpi_type(local_shape, negative_offset, positive_offset, dirs)
+
+                logger.debug(
+                    f"Receiving fractal surface from rank {rank}, MPI type={mpi_type.decode()}"
+                )
+                request = self.comm.Irecv([self.fractalsurface, mpi_type], rank)
+                requests.append(request)
+
+        if len(requests) > 0:
+            requests[0].Waitall(requests)
+
+        # Update start and stop to local bounds
+        self.start = np.maximum(self.start, 0)
+        self.start = np.minimum(self.start, self.upper_bound)
+        self.stop = np.maximum(self.stop, 0)
+        self.stop = np.minimum(self.stop, self.upper_bound)
+
+        logger.debug(
+            f"Generated fractal surface: start={self.start}, stop={self.stop}, size={self.size}, fractalrange={self.fractalrange}"
+        )
+
+    def calculate_starts_and_subshape(
+        self,
+        shape: npt.NDArray[np.int32],
+        negative_offset: npt.NDArray[np.int32],
+        positive_offset: npt.NDArray[np.int32],
+        dirs: npt.NDArray[np.int32],
+        sending: bool = False,
+    ) -> Tuple[npt.NDArray[np.int32], npt.NDArray[np.int32]]:
+        negative_offset = np.where(
+            dirs == Dir.NONE,
+            np.maximum(negative_offset, 0),
+            np.abs(negative_offset),
+        )
+
+        positive_offset = np.where(
+            dirs == Dir.NONE,
+            np.maximum(positive_offset, 0),
+            np.abs(positive_offset),
+        )
+
+        starts = np.select(
+            [dirs == Dir.NEG, dirs == Dir.POS],
+            [0, shape - positive_offset - sending],
+            default=negative_offset,
+        )
+
+        subshape = np.select(
+            [dirs == Dir.NEG, dirs == Dir.POS],
+            [negative_offset + sending, positive_offset + sending],
+            default=shape - negative_offset - positive_offset,
+        )
+
+        return starts, subshape
+
+    def create_mpi_type(
+        self,
+        shape: npt.NDArray[np.int32],
+        negative_offset: npt.NDArray[np.int32],
+        positive_offset: npt.NDArray[np.int32],
+        dirs: npt.NDArray[np.int32],
+        sending: bool = False,
+    ) -> MPI.Datatype:
+        starts, subshape = self.calculate_starts_and_subshape(
+            shape, negative_offset, positive_offset, dirs, sending
+        )
+
+        mpi_type = MPI.FLOAT.Create_subarray(shape.tolist(), subshape.tolist(), starts.tolist())
+        mpi_type.Commit()
+        return mpi_type

gprMax/user_objects/cmds_geometry/add_surface_roughness.py

@@ -20,8 +20,9 @@ import logging
 
 import numpy as np
 
-from gprMax.fractals.fractal_surface import FractalSurface
+from gprMax.fractals.fractal_surface import FractalSurface, MPIFractalSurface
 from gprMax.grid.fdtd_grid import FDTDGrid
+from gprMax.grid.mpi_grid import MPIGrid
 from gprMax.user_objects.rotatable import RotatableMixin
 from gprMax.user_objects.user_objects import GeometryUserObject
 from gprMax.utilities.utilities import round_value
@@ -37,18 +38,19 @@ class AddSurfaceRoughness(RotatableMixin, GeometryUserObject):
     Attributes:
         p1: list of the lower left (x,y,z) coordinates of a surface on a
             FractalBox class.
-        p2: list of the upper right (x,y,z) coordinates of a surface on a
-            FractalBox class.
-        frac_dim: float for the fractal dimension which, for an orthogonal
-                    parallelepiped, should take values between zero and three.
-        weighting: list with weightings in the first and second direction of
-                    the surface.
-        limits: ist to define lower and upper limits for a range over which
-                    the surface roughness can vary.
-        fractal_box_id: string identifier for the FractalBox class
-                        that the surface roughness should be applied to.
-        seed: (optional) float parameter which controls the seeding of the random
-                number generator used to create the fractals.
+        p2: list of the upper right (x,y,z) coordinates of a surface on
+            a FractalBox class.
+        frac_dim: float for the fractal dimension which, for an
+            orthogonal parallelepiped, should take values between zero
+            and three.
+        weighting: list with weightings in the first and second
+            direction of the surface.
+        limits: list to define lower and upper limits for a range over
+            which the surface roughness can vary.
+        fractal_box_id: string identifier for the FractalBox class that
+            the surface roughness should be applied to.
+        seed: (optional) float parameter which controls the seeding of
+            the random number generator used to create the fractals.
     """
 
     @property
@@ -95,145 +97,116 @@ class AddSurfaceRoughness(RotatableMixin, GeometryUserObject):
         if volumes:
             volume = volumes[0]
         else:
-            logger.exception(f"{self.__str__()} cannot find FractalBox {fractal_box_id}")
-            raise ValueError
+            raise ValueError(f"{self.__str__()} cannot find FractalBox {fractal_box_id}")
 
         uip = self._create_uip(grid)
-        _, p1, p2 = uip.check_box_points(p1, p2, self.__str__())
-        xs, ys, zs = p1
-        xf, yf, zf = p2
+        discretised_p1, discretised_p2 = uip.check_output_object_bounds(p1, p2, self.__str__())
+        xs, ys, zs = discretised_p1
+        xf, yf, zf = discretised_p2
 
         if frac_dim < 0:
-            logger.exception(
+            raise ValueError(
                 f"{self.__str__()} requires a positive value for the fractal dimension"
             )
-            raise ValueError
         if weighting[0] < 0:
-            logger.exception(
-                f"{self.__str__()} requires a positive value for the "
-                "fractal weighting in the first direction of the surface"
+            raise ValueError(
+                f"{self.__str__()} requires a positive value for the fractal weighting in the first"
+                " direction of the surface"
             )
-            raise ValueError
         if weighting[1] < 0:
-            logger.exception(
-                f"{self.__str__()} requires a positive value for the "
-                "fractal weighting in the second direction of the surface"
+            raise ValueError(
+                f"{self.__str__()} requires a positive value for the fractal weighting in the"
+                " second direction of the surface"
             )
-            raise ValueError
 
         # Check for valid orientations
+        if np.count_nonzero(discretised_p1 == discretised_p2) != 1:
+            raise ValueError(f"{self.__str__()} dimensions are not specified correctly")
+
         if xs == xf:
-            if ys == yf or zs == zf:
-                logger.exception(f"{self.__str__()} dimensions are not specified correctly")
-                raise ValueError
-            if xs not in [volume.xs, volume.xf]:
-                logger.exception(
-                    f"{self.__str__()} can only be used on the external surfaces of a fractal box"
-                )
-                raise ValueError
-            fractalrange = (
-                round_value(limits[0] / grid.dx),
-                round_value(limits[1] / grid.dx),
-            )
             # xminus surface
             if xs == volume.xs:
-                if fractalrange[0] < 0 or fractalrange[1] > volume.xf:
-                    logger.exception(
-                        f"{self.__str__()} cannot apply fractal surface "
-                        "to fractal box as it would exceed either the "
-                        "upper coordinates of the fractal box or the "
-                        "domain in the x direction"
-                    )
-                    raise ValueError
+                lower_bound = discretised_p1
+                upper_bound = uip.discretise_point((limits[1], p2[1], p2[2]))
+                grid_bound = uip.discretise_point((limits[0], p1[1], p1[2]))
+                fractalrange = (grid_bound[0], upper_bound[0])
                 requestedsurface = "xminus"
             # xplus surface
             elif xf == volume.xf:
-                if fractalrange[0] < volume.xs or fractalrange[1] > grid.nx:
-                    logger.exception(
-                        f"{self.__str__()} cannot apply fractal surface "
-                        "to fractal box as it would exceed either the "
-                        "lower coordinates of the fractal box or the "
-                        "domain in the x direction"
-                    )
-                    raise ValueError
+                lower_bound = uip.discretise_point((limits[0], p1[1], p1[2]))
+                upper_bound = discretised_p2
+                grid_bound = uip.discretise_point((limits[1], p2[1], p2[2]))
+                fractalrange = (lower_bound[0], grid_bound[0])
                 requestedsurface = "xplus"
-
-        elif ys == yf:
-            if zs == zf:
-                logger.exception(f"{self.__str__()} dimensions are not specified correctly")
-                raise ValueError
-            if ys not in [volume.ys, volume.yf]:
-                logger.exception(
-                    f"{self.__str__()} can only be used on the external "
-                    + "surfaces of a fractal box"
+            else:
+                raise ValueError(
+                    f"{self.__str__()} can only be used on the external surfaces of a fractal box"
                 )
-                raise ValueError
-            fractalrange = (
-                round_value(limits[0] / grid.dy),
-                round_value(limits[1] / grid.dy),
-            )
+        elif ys == yf:
             # yminus surface
             if ys == volume.ys:
-                if fractalrange[0] < 0 or fractalrange[1] > volume.yf:
-                    logger.exception(
-                        f"{self.__str__()} cannot apply fractal surface "
-                        "to fractal box as it would exceed either the "
-                        "upper coordinates of the fractal box or the "
-                        "domain in the y direction"
-                    )
-                    raise ValueError
+                lower_bound = discretised_p1
+                upper_bound = uip.discretise_point((p2[0], limits[1], p2[2]))
+                grid_bound = uip.discretise_point((p1[0], limits[0], p1[2]))
+                fractalrange = (grid_bound[1], upper_bound[1])
                 requestedsurface = "yminus"
             # yplus surface
             elif yf == volume.yf:
-                if fractalrange[0] < volume.ys or fractalrange[1] > grid.ny:
-                    logger.exception(
-                        f"{self.__str__()} cannot apply fractal surface "
-                        "to fractal box as it would exceed either the "
-                        "lower coordinates of the fractal box or the "
-                        "domain in the y direction"
-                    )
-                    raise ValueError
+                lower_bound = uip.discretise_point((p1[0], limits[0], p1[2]))
+                upper_bound = discretised_p2
+                grid_bound = uip.discretise_point((p2[0], limits[1], p2[2]))
+                fractalrange = (lower_bound[1], grid_bound[1])
                 requestedsurface = "yplus"
-
-        elif zs == zf:
-            if zs not in [volume.zs, volume.zf]:
-                logger.exception(
-                    f"{self.__str__()} can only be used on the external "
-                    + "surfaces of a fractal box"
+            else:
+                raise ValueError(
+                    f"{self.__str__()} can only be used on the external surfaces of a fractal box"
                 )
-                raise ValueError
-            fractalrange = (
-                round_value(limits[0] / grid.dz),
-                round_value(limits[1] / grid.dz),
-            )
+        elif zs == zf:
             # zminus surface
             if zs == volume.zs:
-                if fractalrange[0] < 0 or fractalrange[1] > volume.zf:
-                    logger.exception(
-                        f"{self.__str__()} cannot apply fractal surface "
-                        "to fractal box as it would exceed either the "
-                        "upper coordinates of the fractal box or the "
-                        "domain in the x direction"
-                    )
-                    raise ValueError
+                lower_bound = discretised_p1
+                upper_bound = uip.discretise_point((p2[0], p2[1], limits[1]))
+                grid_bound = uip.discretise_point((p1[0], p1[1], limits[0]))
+                fractalrange = (grid_bound[2], upper_bound[2])
                 requestedsurface = "zminus"
             # zplus surface
             elif zf == volume.zf:
-                if fractalrange[0] < volume.zs or fractalrange[1] > grid.nz:
-                    logger.exception(
-                        f"{self.__str__()} cannot apply fractal surface "
-                        "to fractal box as it would exceed either the "
-                        "lower coordinates of the fractal box or the "
-                        "domain in the z direction"
-                    )
-                    raise ValueError
+                lower_bound = uip.discretise_point((p1[0], p1[1], limits[0]))
+                upper_bound = discretised_p2
+                grid_bound = uip.discretise_point((p2[0], p2[1], limits[1]))
+                fractalrange = (lower_bound[2], grid_bound[2])
                 requestedsurface = "zplus"
-
+            else:
+                raise ValueError(
+                    f"{self.__str__()} can only be used on the external surfaces of a fractal box"
+                )
         else:
-            logger.exception(f"{self.__str__()} dimensions are not specified correctly")
-            raise ValueError
+            raise ValueError(f"{self.__str__()} dimensions are not specified correctly")
 
-        surface = FractalSurface(xs, xf, ys, yf, zs, zf, frac_dim, seed)
+        if any(lower_bound < volume.start):
+            raise ValueError(
+                f"{self.__str__()} cannot apply fractal surface to"
+                " fractal box as it would exceed the lower coordinates"
+                " of the fractal box."
+            )
+        elif any(upper_bound > volume.stop):
+            raise ValueError(
+                f"{self.__str__()} cannot apply fractal surface to"
+                " fractal box as it would exceed the upper coordinates"
+                " of the fractal box."
+            )
+
+        # Check lower or upper extent of the fractal surface (depending
+        # if the fractal surface has been applied in the minus or plus
+        # direction).
+        uip.point_within_bounds(grid_bound, f"{self.__str__()}")
+
+        if isinstance(grid, MPIGrid):
+            surface = MPIFractalSurface(
+                xs, xf, ys, yf, zs, zf, frac_dim, seed, grid.comm, grid.nx, grid.ny, grid.nz
+            )
+        else:
+            surface = FractalSurface(xs, xf, ys, yf, zs, zf, frac_dim, seed)
         surface.surfaceID = requestedsurface
         surface.fractalrange = fractalrange
         surface.operatingonID = volume.ID
@@ -242,18 +215,21 @@ class AddSurfaceRoughness(RotatableMixin, GeometryUserObject):
         # List of existing surfaces IDs
         existingsurfaceIDs = [x.surfaceID for x in volume.fractalsurfaces]
         if surface.surfaceID in existingsurfaceIDs:
-            logger.exception(
+            raise ValueError(
                 f"{self.__str__()} has already been used on the {surface.surfaceID} surface"
             )
-            raise ValueError
 
         surface.generate_fractal_surface()
+
         volume.fractalsurfaces.append(surface)
 
+        p3 = uip.round_to_grid_static_point(p1)
+        p4 = uip.round_to_grid_static_point(p2)
+
         logger.info(
-            f"{self.grid_name(grid)}Fractal surface from {xs * grid.dx:g}m, "
-            f"{ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, "
-            f"{yf * grid.dy:g}m, {zf * grid.dz:g}m with fractal dimension "
+            f"{self.grid_name(grid)}Fractal surface from {p3[0]:g}m, "
+            f"{p3[1]:g}m, {p3[2]:g}m, to {p4[0]:g}m, "
+            f"{p4[1]:g}m, {p4[2]:g}m with fractal dimension "
             f"{surface.dimension:g}, fractal weightings {surface.weighting[0]:g}, "
             f"{surface.weighting[1]:g}, fractal seeding {surface.seed}, "
             f"and range {limits[0]:g}m to {limits[1]:g}m, added to "

gprMax/user_objects/cmds_geometry/fractal_box.py

@@ -252,7 +252,8 @@ class FractalBox(RotatableMixin, GeometryUserObject):
                     )
                     self.volume.fractalvolume *= materialnumID
                 else:
-                    self.volume.generate_fractal_volume()
+                    if not self.volume.generate_fractal_volume():
+                        return
                     for i in range(0, self.volume.nx):
                         for j in range(0, self.volume.ny):
                             for k in range(0, self.volume.nz):
@@ -268,6 +269,10 @@ class FractalBox(RotatableMixin, GeometryUserObject):
                 # TODO: Allow extract of rough surface profile (to print/file?)
                 for surface in self.volume.fractalsurfaces:
                     if surface.surfaceID == "xminus":
+                        surface.fractalrange = (
+                            max(surface.fractalrange[0], 0),
+                            min(surface.fractalrange[1], grid.nx),
+                        )
                         for i in range(surface.fractalrange[0], surface.fractalrange[1]):
                             for j in range(surface.ys, surface.yf):
                                 for k in range(surface.zs, surface.zf):
@@ -291,6 +296,10 @@ class FractalBox(RotatableMixin, GeometryUserObject):
                                         ] = 0
 
                     elif surface.surfaceID == "xplus":
+                        surface.fractalrange = (
+                            max(surface.fractalrange[0], 0),
+                            min(surface.fractalrange[1], grid.nx),
+                        )
                         if not surface.ID:
                             for i in range(surface.fractalrange[0], surface.fractalrange[1]):
                                 for j in range(surface.ys, surface.yf):
@@ -399,6 +408,10 @@ class FractalBox(RotatableMixin, GeometryUserObject):
                                         root += 1
 
                     elif surface.surfaceID == "yminus":
+                        surface.fractalrange = (
+                            max(surface.fractalrange[0], 0),
+                            min(surface.fractalrange[1], grid.ny),
+                        )
                         for i in range(surface.xs, surface.xf):
                             for j in range(surface.fractalrange[0], surface.fractalrange[1]):
                                 for k in range(surface.zs, surface.zf):
@@ -422,6 +435,10 @@ class FractalBox(RotatableMixin, GeometryUserObject):
                                         ] = 0
 
                     elif surface.surfaceID == "yplus":
+                        surface.fractalrange = (
+                            max(surface.fractalrange[0], 0),
+                            min(surface.fractalrange[1], grid.ny),
+                        )
                         if not surface.ID:
                             for i in range(surface.xs, surface.xf):
                                 for j in range(surface.fractalrange[0], surface.fractalrange[1]):
@@ -530,6 +547,10 @@ class FractalBox(RotatableMixin, GeometryUserObject):
                                         root += 1
 
                     elif surface.surfaceID == "zminus":
+                        surface.fractalrange = (
+                            max(surface.fractalrange[0], 0),
+                            min(surface.fractalrange[1], grid.nz),
+                        )
                         for i in range(surface.xs, surface.xf):
                             for j in range(surface.ys, surface.yf):
                                 for k in range(surface.fractalrange[0], surface.fractalrange[1]):
@@ -553,6 +574,10 @@ class FractalBox(RotatableMixin, GeometryUserObject):
                                         ] = 0
 
                     elif surface.surfaceID == "zplus":
+                        surface.fractalrange = (
+                            max(surface.fractalrange[0], 0),
+                            min(surface.fractalrange[1], grid.nz),
+                        )
                         if not surface.ID:
                             for i in range(surface.xs, surface.xf):
                                 for j in range(surface.ys, surface.yf):