From 14a4c342ac3ad27c488bac2b8890c94885f26bf1 Mon Sep 17 00:00:00 2001 From: craig-warren Date: Mon, 21 Oct 2019 20:46:03 +0100 Subject: [PATCH] Working on getting heterogeneous_soil example to run, i.e. with fractals & dispersive materials. --- gprMax/fractals.py | 8 +++--- gprMax/grid.py | 55 +++++++++++++++++++++++++++------------ gprMax/materials.py | 12 ++++----- gprMax/model_build_run.py | 9 +++++-- gprMax/pml.py | 2 +- 5 files changed, 57 insertions(+), 29 deletions(-) diff --git a/gprMax/fractals.py b/gprMax/fractals.py index 7352e056..24b3d836 100644 --- a/gprMax/fractals.py +++ b/gprMax/fractals.py @@ -77,7 +77,7 @@ class FractalSurface: elif self.zs == self.zf: surfacedims = (self.nx, self.ny) - self.fractalsurface = np.zeros(surfacedims, dtype=config.dtypes['complex']) + self.fractalsurface = np.zeros(surfacedims, dtype=config.sim_config.dtypes['complex']) # Positional vector at centre of array, scaled by weighting v1 = np.array([self.weighting[0] * (surfacedims[0]) / 2, self.weighting[1] @@ -93,7 +93,7 @@ class FractalSurface: A = fftpack.fftshift(A) # Generate fractal - generate_fractal2D(surfacedims[0], surfacedims[1], config.hostinfo['ompthreads'], + generate_fractal2D(surfacedims[0], surfacedims[1], config.model_configs[G.model_num].ompthreads, self.b, self.weighting, v1, A, self.fractalsurface) # Shift the zero frequency component to start of the array self.fractalsurface = fftpack.ifftshift(self.fractalsurface) @@ -164,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.dtypes['complex']) + self.fractalvolume = np.zeros((self.nx, self.ny, self.nz), dtype=config.sim_config.dtypes['complex']) # Positional vector at centre of array, scaled by weighting v1 = np.array([self.weighting[0] * self.nx / 2, self.weighting[1] @@ -180,7 +180,7 @@ class FractalVolume: A = fftpack.fftshift(A) # Generate fractal - generate_fractal3D(self.nx, self.ny, self.nz, config.hostinfo['ompthreads'], + generate_fractal3D(self.nx, self.ny, self.nz, config.model_configs[G.model_num].ompthreads, self.b, self.weighting, v1, A, self.fractalvolume) # Shift the zero frequency component to the start of the array diff --git a/gprMax/grid.py b/gprMax/grid.py index d8cc4b2d..e880f9ab 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -148,12 +148,18 @@ class FDTDGrid: def initialise_field_arrays(self): """Initialise arrays for the electric and magnetic field components.""" - self.Ex = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes['float_or_double']) - self.Ey = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes['float_or_double']) - self.Ez = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes['float_or_double']) - self.Hx = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes['float_or_double']) - self.Hy = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes['float_or_double']) - self.Hz = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.sim_config.dtypes['float_or_double']) + self.Ex = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), + dtype=config.sim_config.dtypes['float_or_double']) + self.Ey = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), + dtype=config.sim_config.dtypes['float_or_double']) + self.Ez = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), + dtype=config.sim_config.dtypes['float_or_double']) + self.Hx = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), + dtype=config.sim_config.dtypes['float_or_double']) + self.Hy = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), + dtype=config.sim_config.dtypes['float_or_double']) + self.Hz = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), + dtype=config.sim_config.dtypes['float_or_double']) def initialise_grids(self): """Initialise all grids.""" @@ -163,15 +169,26 @@ class FDTDGrid: def initialise_std_update_coeff_arrays(self): """Initialise arrays for storing update coefficients.""" - self.updatecoeffsE = np.zeros((len(self.materials), 5), dtype=config.sim_config.dtypes['float_or_double']) - self.updatecoeffsH = np.zeros((len(self.materials), 5), dtype=config.sim_config.dtypes['float_or_double']) + self.updatecoeffsE = np.zeros((len(self.materials), 5), + dtype=config.sim_config.dtypes['float_or_double']) + self.updatecoeffsH = np.zeros((len(self.materials), 5), + dtype=config.sim_config.dtypes['float_or_double']) def initialise_dispersive_arrays(self, dtype): - """Initialise arrays for storing coefficients when there are dispersive materials present.""" - self.Tx = np.zeros((config.materials['maxpoles'], self.nx + 1, self.ny + 1, self.nz + 1), dtype=dtype) - self.Ty = np.zeros((config.materials['maxpoles'], self.nx + 1, self.ny + 1, self.nz + 1), dtype=dtype) - self.Tz = np.zeros((config.materials['maxpoles'], self.nx + 1, self.ny + 1, self.nz + 1), dtype=dtype) - self.updatecoeffsdispersive = np.zeros((len(self.materials), 3 * config.model_configs[self.model_num].materials['maxpoles']), dtype=dtype) + """Initialise arrays for storing coefficients when there are dispersive materials present. + + Args: + dtype (dtype): Dtype to use for dispersive arrays. + """ + self.Tx = np.zeros((config.model_configs[self.model_num].materials['maxpoles'], + self.nx + 1, self.ny + 1, self.nz + 1), dtype=dtype) + self.Ty = np.zeros((config.model_configs[self.model_num].materials['maxpoles'], + self.nx + 1, self.ny + 1, self.nz + 1), dtype=dtype) + self.Tz = np.zeros((config.model_configs[self.model_num].materials['maxpoles'], + self.nx + 1, self.ny + 1, self.nz + 1), dtype=dtype) + self.updatecoeffsdispersive = np.zeros((len(self.materials), 3 * + config.model_configs[self.model_num].materials['maxpoles']), + dtype=dtype) def mem_est_basic(self): """Estimate the amount of memory (RAM) required for grid arrays. @@ -220,7 +237,7 @@ class FDTDGrid: """ mem_use = int(3 * config.model_configs[self.model_num].materials['maxpoles'] * - (G.nx + 1) * (G.ny + 1) * (G.nz + 1) * + (self.nx + 1) * (self.ny + 1) * (self.nz + 1) * np.dtype(config.model_configs[self.model_num].materials['dispersivedtype']).itemsize) return mem_use @@ -309,11 +326,17 @@ class CUDAGrid(FDTDGrid): self.Hy_gpu = gpuarray.to_gpu(self.Hy) self.Hz_gpu = gpuarray.to_gpu(self.Hz) - def initialise_dispersive_arrays(self): - """Initialise dispersive material coefficient arrays on GPU.""" + def initialise_dispersive_arrays(self, dtype): + """Initialise dispersive material coefficient arrays on GPU. + + Args: + dtype (dtype): Dtype to use for dispersive arrays. + """ import pycuda.gpuarray as gpuarray + super().initialise_dispersive_arrays(dtype) + self.Tx_gpu = gpuarray.to_gpu(self.Tx) self.Ty_gpu = gpuarray.to_gpu(self.Ty) self.Tz_gpu = gpuarray.to_gpu(self.Tz) diff --git a/gprMax/materials.py b/gprMax/materials.py index fe5689f6..b12d35e0 100644 --- a/gprMax/materials.py +++ b/gprMax/materials.py @@ -132,17 +132,17 @@ class DispersiveMaterial(Material): # The implementation of the dispersive material modelling comes from the # derivation in: http://dx.doi.org/10.1109/TAP.2014.2308549 self.w = np.zeros(config.model_configs[G.model_num].materials['maxpoles'], - dtype=config.materials['dispersivedtype']) + dtype=config.model_configs[G.model_num].materials['dispersivedtype']) self.q = np.zeros(config.model_configs[G.model_num].materials['maxpoles'], - dtype=config.materials['dispersivedtype']) + dtype=config.model_configs[G.model_num].materials['dispersivedtype']) self.zt = np.zeros(config.model_configs[G.model_num].materials['maxpoles'], - dtype=config.materials['dispersivedtype']) + dtype=config.model_configs[G.model_num].materials['dispersivedtype']) self.zt2 = np.zeros(config.model_configs[G.model_num].materials['maxpoles'], - dtype=config.materials['dispersivedtype']) + dtype=config.model_configs[G.model_num].materials['dispersivedtype']) self.eqt = np.zeros(config.model_configs[G.model_num].materials['maxpoles'], - dtype=config.materials['dispersivedtype']) + dtype=config.model_configs[G.model_num].materials['dispersivedtype']) self.eqt2 = np.zeros(config.model_configs[G.model_num].materials['maxpoles'], - dtype=config.materials['dispersivedtype']) + dtype=config.model_configs[G.model_num].materials['dispersivedtype']) for x in range(self.poles): if 'debye' in self.type: diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py index 5bd08364..b9a37089 100644 --- a/gprMax/model_build_run.py +++ b/gprMax/model_build_run.py @@ -86,6 +86,12 @@ class ModelBuildRun: # Monitor memory usage self.p = None + # Set number of OpenMP threads to physical threads at this point to be + # used with threaded model building methods, e.g. fractals. Can be + # changed by #num_threads command in input file or via API later for + # use with CPU solver. + config.model_configs[self.G.model_num].ompthreads = set_omp_threads(config.model_configs[self.G.model_num].ompthreads) + def build(self): """Builds the Yee cells for a model.""" @@ -305,9 +311,8 @@ class ModelBuildRun: self.create_output_directory() log.info(f'\nOutput file: {config.model_configs[self.G.model_num].output_file_path_ext}') - # Set and check number of OpenMP threads + # Check number of OpenMP threads if config.sim_config.general['cpu']: - config.model_configs[self.G.model_num].ompthreads = set_omp_threads(config.model_configs[self.G.model_num].ompthreads) log.info(f'CPU (OpenMP) threads for solving: {config.model_configs[self.G.model_num].ompthreads}\n') if config.model_configs[self.G.model_num].ompthreads > config.sim_config.hostinfo['physicalcores']: log.warning(Fore.RED + f"You have specified more threads ({config.model_configs[self.G.model_num].ompthreads}) than available physical CPU cores ({config.sim_config.hostinfo['physicalcores']}). This may lead to degraded performance." + Style.RESET_ALL) diff --git a/gprMax/pml.py b/gprMax/pml.py index f5ad1bf1..08cb8408 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.py @@ -238,7 +238,7 @@ class PML: self.HPhi1 = np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz), dtype=config.sim_config.dtypes['float_or_double']) self.HPhi2 = np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz), - dtype=config.sim_config.sim_config.dtypes['float_or_double']) + dtype=config.sim_config.dtypes['float_or_double']) def calculate_update_coeffs(self, er, mr, G): """Calculates electric and magnetic update coefficients for the PML.