Working on getting heterogeneous_soil example to run, i.e. with fractals & dispersive materials.

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

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)

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:

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)

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.