From 12589aae272c594ac36575de7594c10bf8429d4a Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Wed, 20 Apr 2022 10:23:30 +0100 Subject: [PATCH] Correct fractal float/double issue, and update random number generation for grass. --- gprMax/cmds_geometry/add_grass.py | 2 +- gprMax/cmds_geometry/add_surface_roughness.py | 2 +- gprMax/cmds_geometry/fractal_box_builder.py | 4 +- gprMax/cython/fractals_generate.pyx | 24 ++++---- gprMax/fractals.py | 61 ++++++++++--------- 5 files changed, 46 insertions(+), 47 deletions(-) diff --git a/gprMax/cmds_geometry/add_grass.py b/gprMax/cmds_geometry/add_grass.py index 1d7ab457..762e4b61 100644 --- a/gprMax/cmds_geometry/add_grass.py +++ b/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 diff --git a/gprMax/cmds_geometry/add_surface_roughness.py b/gprMax/cmds_geometry/add_surface_roughness.py index d33d61c1..47db0a13 100644 --- a/gprMax/cmds_geometry/add_surface_roughness.py +++ b/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}.') diff --git a/gprMax/cmds_geometry/fractal_box_builder.py b/gprMax/cmds_geometry/fractal_box_builder.py index 5d88cacd..35946b9e 100644 --- a/gprMax/cmds_geometry/fractal_box_builder.py +++ b/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') diff --git a/gprMax/cython/fractals_generate.pyx b/gprMax/cython/fractals_generate.pyx index d300f875..0c59622a 100644 --- a/gprMax/cython/fractals_generate.pyx +++ b/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) - - # Catch potential divide by zero - if rr == 0: - rr = 0.9 - - fractalvolume[i, j, k] = A[i, j, k] * 1 / (rr**b) + B = rr**b + if B == 0: + B = 0.9 + + fractalvolume[i, j, k] = A[i, j, k] / B diff --git a/gprMax/fractals.py b/gprMax/fractals.py index d908b7ad..19a8a7e3 100644 --- a/gprMax/fractals.py +++ 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])