Correct fractal float/double issue, and update random number generation for grass.

gprMax/cmds_geometry/add_grass.py

@@ -173,7 +173,7 @@ class AddGrass(UserObjectGeometry):
         # Set the fractal range to scale the fractal distribution between zero and one
         surface.fractalrange = (0, 1)
         surface.operatingonID = volume.ID
-        surface.generate_fractal_surface(grid)
+        surface.generate_fractal_surface()
         if n_blades > surface.fractalsurface.shape[0] * surface.fractalsurface.shape[1]:
             logger.exception(self.__str__() + ' the specified surface is not large enough for the number of grass blades/roots specified')
             raise ValueError

gprMax/cmds_geometry/add_surface_roughness.py

@@ -191,7 +191,7 @@ class AddSurfaceRoughness(UserObjectGeometry):
             logger.exception(self.__str__() + f' has already been used on the {surface.surfaceID} surface')
             raise ValueError
 
-        surface.generate_fractal_surface(grid)
+        surface.generate_fractal_surface()
         volume.fractalsurfaces.append(surface)
 
         logger.info(self.grid_name(grid) + f'Fractal surface from {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m with fractal dimension {surface.dimension:g}, fractal weightings {surface.weighting[0]:g}, {surface.weighting[1]:g}, fractal seeding {surface.seed}, and range {limits[0]:g}m to {limits[1]:g}m, added to {surface.operatingonID}.')

gprMax/cmds_geometry/fractal_box_builder.py

@@ -80,7 +80,7 @@ class FractalBoxBuilder(UserObjectGeometry):
                     materialnumID = next(x.numID for x in grid.materials if x.ID == volume.operatingonID)
                     volume.fractalvolume *= materialnumID
                 else:
-                    volume.generate_fractal_volume(grid)
+                    volume.generate_fractal_volume()
                     volume.fractalvolume += volume.mixingmodel.startmaterialnum
 
                 volume.generate_volume_mask()
@@ -306,7 +306,7 @@ class FractalBoxBuilder(UserObjectGeometry):
                     logger.exception(self.__str__() + ' is being used with a single material and no modifications, therefore please use a #box command instead.')
                     raise ValueError
                 else:
-                    volume.generate_fractal_volume(grid)
+                    volume.generate_fractal_volume()
                     volume.fractalvolume += volume.mixingmodel.startmaterialnum
 
                 data = volume.fractalvolume.astype('int16', order='C')

gprMax/cython/fractals_generate.pyx

@@ -37,7 +37,7 @@ cpdef void generate_fractal2D(int nx, int ny, int nthreads, int b, np.float64_t[
     """
 
     cdef Py_ssize_t i, j
-    cdef float v2x, v2y, rr
+    cdef double v2x, v2y, rr, B
 
     for i in prange(nx, nogil=True, schedule='static', num_threads=nthreads):
         for j in range(ny):
@@ -47,12 +47,11 @@ cpdef void generate_fractal2D(int nx, int ny, int nthreads, int b, np.float64_t[
 
                 # Calulate norm of v2 - v1
                 rr = ((v2x - v1[0])**2 + (v2y - v1[1])**2)**(1/2)
+                B = rr**b
+                if B == 0:
+                    B = 0.9
 
-                # Catch potential divide by zero
-                if rr == 0:
-                    rr = 0.9
-
-                fractalsurface[i, j] = A[i, j] * 1 / (rr**b)
+                fractalsurface[i, j] = A[i, j] / B
 
 
 cpdef void generate_fractal3D(int nx, int ny, int nz, int nthreads, int b, np.float64_t[:] weighting, np.float64_t[:] v1, np.complex128_t[:, :, ::1] A, float_or_double_complex[:, :, ::1] fractalvolume):
@@ -69,7 +68,7 @@ cpdef void generate_fractal3D(int nx, int ny, int nz, int nthreads, int b, np.fl
     """
 
     cdef Py_ssize_t i, j, k
-    cdef float v2x, v2y, v2z, rr
+    cdef double v2x, v2y, v2z, rr, B
 
     for i in prange(nx, nogil=True, schedule='static', num_threads=nthreads):
         for j in range(ny):
@@ -81,9 +80,8 @@ cpdef void generate_fractal3D(int nx, int ny, int nz, int nthreads, int b, np.fl
 
                 # Calulate norm of v2 - v1
                 rr = ((v2x - v1[0])**2 + (v2y - v1[1])**2 + (v2z - v1[2])**2)**(1/2)
+                B = rr**b
+                if B == 0:
+                    B = 0.9 
                     
-                # Catch potential divide by zero
-                if rr == 0:
-                    rr = 0.9
-
-                fractalvolume[i, j, k] = A[i, j, k] * 1 / (rr**b)
+                fractalvolume[i, j, k] = A[i, j, k] / B

gprMax/fractals.py

@@ -61,13 +61,8 @@ class FractalSurface:
         self.filldepth = 0
         self.grass = []
 
-    def generate_fractal_surface(self, G):
-        """Generate a 2D array with a fractal distribution.
-
-        Args:
-            G (class): Grid class instance - holds essential parameters
-                        describing the model.
-        """
+    def generate_fractal_surface(self):
+        """Generate a 2D array with a fractal distribution."""
 
         if self.xs == self.xf:
             surfacedims = (self.ny, self.nz)
@@ -76,7 +71,7 @@ class FractalSurface:
         elif self.zs == self.zf:
             surfacedims = (self.nx, self.ny)
 
-        self.fractalsurface = np.zeros(surfacedims, dtype=config.sim_config.dtypes['complex'])
+        self.fractalsurface = np.zeros(surfacedims, dtype=np.complex128)
 
         # Positional vector at centre of array, scaled by weighting
         v1 = np.array([self.weighting[0] * (surfacedims[0]) / 2, self.weighting[1]
@@ -96,8 +91,13 @@ class FractalSurface:
                            self.b, self.weighting, v1, A, self.fractalsurface)
         # Shift the zero frequency component to start of the array
         self.fractalsurface = fftpack.ifftshift(self.fractalsurface)
-        # Take the real part (numerical errors can give rise to an imaginary part) of the IFFT
-        self.fractalsurface = np.real(fftpack.ifftn(self.fractalsurface))
+        # Set DC component of FFT to zero
+        self.fractalsurface[0, 0] = 0
+        # 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
+        self.fractalsurface = np.real(fftpack.ifftn(self.fractalsurface)).astype(config.sim_config.dtypes['float_or_double'], 
+                                                                                 copy=False)
         # Scale the fractal volume according to requested range
         fractalmin = np.amin(self.fractalsurface)
         fractalmax = np.amax(self.fractalsurface)
@@ -145,13 +145,8 @@ class FractalVolume:
         self.nbins = 0
         self.fractalsurfaces = []
 
-    def generate_fractal_volume(self, G):
-        """Generate a 3D volume with a fractal distribution.
-
-        Args:
-            G (class): Grid class instance - holds essential parameters
-                        describing the model.
-        """
+    def generate_fractal_volume(self):
+        """Generate a 3D volume with a fractal distribution."""
 
         # Scale filter according to size of fractal volume
         if self.nx == 1:
@@ -169,7 +164,7 @@ class FractalVolume:
         # Adjust weighting to account for filter scaling
         self.weighting = np.multiply(self.weighting, filterscaling)
 
-        self.fractalvolume = np.zeros((self.nx, self.ny, self.nz), dtype=config.sim_config.dtypes['complex'])
+        self.fractalvolume = np.zeros((self.nx, self.ny, self.nz), dtype=np.complex128)
 
         # Positional vector at centre of array, scaled by weighting
         v1 = np.array([self.weighting[0] * self.nx / 2, self.weighting[1]
@@ -190,8 +185,14 @@ class FractalVolume:
 
         # Shift the zero frequency component to the start of the array
         self.fractalvolume = fftpack.ifftshift(self.fractalvolume)
-        # Take the real part (numerical errors can give rise to an imaginary part) of the IFFT
-        self.fractalvolume = np.real(fftpack.ifftn(self.fractalvolume))
+        # Set DC component of FFT to zero
+        self.fractalvolume[0, 0, 0] = 0
+        # 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
+        self.fractalvolume = np.real(fftpack.ifftn(self.fractalvolume)).astype(config.sim_config.dtypes['float_or_double'], 
+                                                                               copy=False)
+
         # Bin fractal values
         bins = np.linspace(np.amin(self.fractalvolume), np.amax(self.fractalvolume), self.nbins)
         for j in range(self.ny):
@@ -229,16 +230,16 @@ class Grass:
         self.seed = None
 
         # Randomly defined parameters that will be used to calculate geometry
-        self.R1 = np.random.RandomState(self.seed)
-        self.R2 = np.random.RandomState(self.seed)
-        self.R3 = np.random.RandomState(self.seed)
-        self.R4 = np.random.RandomState(self.seed)
-        self.R5 = np.random.RandomState(self.seed)
-        self.R6 = np.random.RandomState(self.seed)
+        self.R1 = np.random.default_rng(seed=self.seed)
+        self.R2 = np.random.default_rng(seed=self.seed)
+        self.R3 = np.random.default_rng(seed=self.seed)
+        self.R4 = np.random.default_rng(seed=self.seed)
+        self.R5 = np.random.default_rng(seed=self.seed)
+        self.R6 = np.random.default_rng(seed=self.seed)
 
         for i in range(self.numblades):
-            self.geometryparams[i, 0] = 10 + 20 * self.R1.random_sample()
-            self.geometryparams[i, 1] = 10 + 20 * self.R2.random_sample()
+            self.geometryparams[i, 0] = 10 + 20 * self.R1.random()
+            self.geometryparams[i, 1] = 10 + 20 * self.R2.random()
             self.geometryparams[i, 2] = self.R3.choice([-1, 1])
             self.geometryparams[i, 3] = self.R4.choice([-1, 1])
 
@@ -273,8 +274,8 @@ class Grass:
             x, y (float): x and y coordinates of grass root.
         """
 
-        self.geometryparams[root, 4] += -1 + 2 * self.R5.random_sample()
-        self.geometryparams[root, 5] += -1 + 2 * self.R6.random_sample()
+        self.geometryparams[root, 4] += -1 + 2 * self.R5.random()
+        self.geometryparams[root, 5] += -1 + 2 * self.R6.random()
         x = round(self.geometryparams[root, 4])
         y = round(self.geometryparams[root, 5])