Code formatting cleanups.

gprMax/config.py

@@ -29,7 +29,7 @@ from .utilities import get_host_info
 # Impedance of free space (Ohms)
 z0 = np.sqrt(m0 / e0)
 
-# General setting for the simulation
+# General setting for the simulation
 #   inputfilepath: path to inputfile location
 #   outputfilepath: path to outputfile location
 #   messages: whether to print all messages as output to stdout or not
@@ -37,7 +37,7 @@ z0 = np.sqrt(m0 / e0)
 #   mode: 2D TMx, 2D TMy, 2D TMz, or 3D
 #   cpu, cuda, opencl: solver type
 general = {'inputfilepath': '', 'outputfilepath': '', 'messages': True,
-            'progressbars': True, 'mode': None, 'cpu': True, 'cuda': False, 'opencl': False}
+           'progressbars': True, 'mode': None, 'cpu': True, 'cuda': False, 'opencl': False}
 
 # Store information about host machine
 hostinfo = get_host_info()
@@ -71,9 +71,9 @@ precision = 'single'
 
 if precision == 'single':
     dtypes = {'float_or_double': np.float32, 'complex': np.complex64,
-            'cython_float_or_double': cython.float, 'cython_complex': cython.floatcomplex,
+              'cython_float_or_double': cython.float, 'cython_complex': cython.floatcomplex,
-            'C_float_or_double': 'float', 'C_complex': 'pycuda::complex<float>'}
+              'C_float_or_double': 'float', 'C_complex': 'pycuda::complex<float>'}
 elif precision == 'double':
     dtypes = {'float_or_double': np.float64, 'complex': np.complex128,
-            'cython_float_or_double': cython.double, 'cython_complex': cython.doublecomplex,
+              'cython_float_or_double': cython.double, 'cython_complex': cython.doublecomplex,
-            'C_float_or_double': 'double', 'C_complex': 'pycuda::complex<double>'}
+              'C_float_or_double': 'double', 'C_complex': 'pycuda::complex<double>'}

gprMax/fields_outputs.py

@@ -44,7 +44,7 @@ def store_outputs(iteration, Ex, Ey, Ez, Hx, Hy, Hz, G):
             else:
                 func = globals()[output]
                 rx.outputs[output][iteration] = func(rx.xcoord, rx.ycoord, rx.zcoord,
-                                                    Hx, Hy, Hz, G)
+                                                     Hx, Hy, Hz, G)
 
     for tl in G.transmissionlines:
         tl.Vtotal[iteration] = tl.voltage[tl.antpos]

gprMax/fractals.py

@@ -81,7 +81,7 @@ class FractalSurface(object):
 
         # Positional vector at centre of array, scaled by weighting
         v1 = np.array([self.weighting[0] * (surfacedims[0]) / 2, self.weighting[1]
-                        * (surfacedims[1]) / 2])
+                       * (surfacedims[1]) / 2])
 
         # 2D array of random numbers to be convolved with the fractal function
         R = np.random.RandomState(self.seed)
@@ -94,7 +94,7 @@ class FractalSurface(object):
 
         # Generate fractal
         generate_fractal2D(surfacedims[0], surfacedims[1], config.hostinfo['ompthreads'],
-                            self.b, self.weighting, v1, A, self.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)
@@ -105,9 +105,9 @@ class FractalSurface(object):
         fractalmax = np.amax(self.fractalsurface)
         fractalrange = fractalmax - fractalmin
         self.fractalsurface = (self.fractalsurface * ((self.fractalrange[1]
-                                - self.fractalrange[0]) / fractalrange)
+                                                       - self.fractalrange[0]) / fractalrange)
-                                + self.fractalrange[0] - ((self.fractalrange[1]
+                               + self.fractalrange[0] - ((self.fractalrange[1]
-                                - self.fractalrange[0]) / fractalrange) * fractalmin)
+                                                          - self.fractalrange[0]) / fractalrange) * fractalmin)
 
 
 class FractalVolume(object):
@@ -169,7 +169,7 @@ class FractalVolume(object):
 
         # Positional vector at centre of array, scaled by weighting
         v1 = np.array([self.weighting[0] * self.nx / 2, self.weighting[1]
-                        * self.ny / 2, self.weighting[2] * self.nz / 2])
+                       * self.ny / 2, self.weighting[2] * self.nz / 2])
 
         # 3D array of random numbers to be convolved with the fractal function
         R = np.random.RandomState(self.seed)
@@ -182,7 +182,7 @@ class FractalVolume(object):
 
         # Generate fractal
         generate_fractal3D(self.nx, self.ny, self.nz, config.hostinfo['ompthreads'],
-                            self.b, self.weighting, v1, A, self.fractalvolume)
+                           self.b, self.weighting, v1, A, self.fractalvolume)
 
         # Shift the zero frequency component to the start of the array
         self.fractalvolume = fftpack.ifftshift(self.fractalvolume)
@@ -221,7 +221,7 @@ class Grass(object):
 
         self.numblades = numblades
         self.geometryparams = np.zeros((self.numblades, 6),
-                                dtype=config.dtypes['float_or_double'])
+                                       dtype=config.dtypes['float_or_double'])
         self.seed = None
 
         # Randomly defined parameters that will be used to calculate geometry
@@ -250,9 +250,9 @@ class Grass(object):
         """
 
         x = (self.geometryparams[blade, 2] * (height / self.geometryparams[blade, 0])
-                * (height / self.geometryparams[blade, 0]))
+             * (height / self.geometryparams[blade, 0]))
         y = (self.geometryparams[blade, 3] * (height / self.geometryparams[blade, 1])
-                * (height / self.geometryparams[blade, 1]))
+             * (height / self.geometryparams[blade, 1]))
         x = round_value(x)
         y = round_value(y)
 

gprMax/gprMax.py

@@ -40,6 +40,7 @@ from .utilities import logo
 from .utilities import open_path_file
 from .utilities import timer
 
+
 def main():
     """This is the main function for gprMax."""
 
@@ -358,8 +359,6 @@ def run_mpi_sim(args, inputfile, usernamespace):
             workobj = dict()
             workobj['currentmodelrun'] = model
             workobj['mpicommname'] = comm.name
-            if optparams:
-                workobj['optparams'] = optparams
             worklist.append(workobj)
         # Add stop sentinels
         worklist += ([StopIteration] * numworkers)

gprMax/input_cmds_file.py

@@ -20,6 +20,7 @@ import os
 import sys
 from io import StringIO
 
+import gprMax.config as config
 from .exceptions import CmdInputError
 
 
@@ -162,7 +163,7 @@ def write_processed_file(processedlines, appendmodelnumber):
         appendmodelnumber (str): Text to append to filename.
     """
 
-    processedfile = os.path.join(outputfilepath, os.path.splitext(inputfilepath)[0] + appendmodelnumber + '_processed.in')
+    processedfile = os.path.join(config.general['outputfilepath'], os.path.splitext(config.general['inputfilepath'])[0] + appendmodelnumber + '_processed.in')
 
     with open(processedfile, 'w') as f:
         for item in processedlines:

gprMax/input_cmds_multiuse.py

@@ -316,7 +316,7 @@ def process_multicmds(multicmds, G):
                 raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires at least six parameters')
 
             # Warn about using a transmission line on GPU
-            if gpus is not None:
+            if config.cuda['gpus'] is not None:
                 raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' A #transmission_line cannot currently be used with GPU solving. Consider using a #voltage_source instead.')
 
             # Check polarity & position parameters

gprMax/materials.py

@@ -196,7 +196,7 @@ class DispersiveMaterial(Material):
         elif 'lorentz' in self.type:
             for pole in range(self.poles):
                 er += ((self.deltaer[pole] * self.tau[pole]**2)
-                        / (self.tau[pole]**2 + 2j * w * self.alpha[pole] - w**2))
+                       / (self.tau[pole]**2 + 2j * w * self.alpha[pole] - w**2))
         elif 'drude' in self.type:
             ersum = 0
             for pole in range(self.poles):
@@ -221,11 +221,11 @@ def process_materials(G):
 
     if config.materials['maxpoles'] == 0:
         materialsdata = [['\nID', '\nName', '\nType', '\neps_r', 'sigma\n[S/m]',
-                            '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']]
+                          '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']]
     else:
         materialsdata = [['\nID', '\nName', '\nType', '\neps_r', 'sigma\n[S/m]',
-                            'Delta\neps_r', 'tau\n[s]', 'omega\n[Hz]', 'delta\n[Hz]',
+                          'Delta\neps_r', 'tau\n[s]', 'omega\n[Hz]', 'delta\n[Hz]',
-                            'gamma\n[Hz]', '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']]
+                          'gamma\n[Hz]', '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']]
 
     for material in G.materials:
         # Calculate update coefficients for specific material
@@ -241,7 +241,7 @@ def process_materials(G):
             z = 0
             for pole in range(config.materials['maxpoles']):
                 G.updatecoeffsdispersive[material.numID, z:z + 3] = (config.e0
-                    * material.eqt2[pole], material.eqt[pole], material.zt[pole])
+                                                                     * material.eqt2[pole], material.eqt[pole], material.zt[pole])
                 z += 3
 
         # Construct information on material properties for printing table
@@ -318,7 +318,7 @@ class PeplinskiSoil(object):
         f = 1.3e9
         w = 2 * np.pi * f
         erealw = DispersiveMaterial.watereri + ((DispersiveMaterial.waterdeltaer)
-                    / (1 + (w * DispersiveMaterial.watertau)**2))
+                                                / (1 + (w * DispersiveMaterial.watertau)**2))
 
         a = 0.65  # Experimentally derived constant
         es = (1.01 + 0.44 * self.rs)**2 - 0.062  #  Relative permittivity of sand particles
@@ -340,7 +340,7 @@ class PeplinskiSoil(object):
         while not muiter.finished:
             # Real part for frequencies in the range 1.4GHz to 18GHz
             er = (1 + (self.rb / self.rs) * ((es**a) - 1) + (muiter[0]**b1 * erealw**a)
-                    - muiter[0]) ** (1 / a)
+                  - muiter[0]) ** (1 / a)
             # Real part for frequencies in the range 0.3GHz to 1.3GHz (linear correction to 1.4-18GHz value)
             er = 1.15 * er - 0.68
 

gprMax/model_build_run.py

@@ -33,6 +33,7 @@ from tqdm import tqdm
 
 import gprMax.config as config
 from .cuda.fields_updates import kernel_template_fields
+from .cuda.snapshots import kernel_template_store_snapshot
 from .cuda.source_updates import kernel_template_sources
 from .cython.fields_updates_normal import update_electric
 from .cython.fields_updates_normal import update_magnetic
@@ -495,19 +496,19 @@ def solve_gpu(currentmodelrun, modelend, G):
 
     # Electric and magnetic field updates - prepare kernels, get kernel functions, and initialise arrays on GPU
     if config.materials['maxpoles'] == 0:
-        kernels_fields = SourceModule(kernels_template_fields.substitute(REAL=config.dtypes['C_float_or_double'], REAL_OR_COMPLEX=config.dtypes['C_complex'], N_updatecoeffsE=G.updatecoeffsE.size, N_updatecoeffsH=G.updatecoeffsH.size, NY_MATCOEFFS=G.updatecoeffsE.shape[1], NY_MATDISPCOEFFS=1, NX_FIELDS=G.Ex.shape[0], NY_FIELDS=G.Ex.shape[1], NZ_FIELDS=G.Ex.shape[2], NX_ID=G.ID.shape[1], NY_ID=G.ID.shape[2], NZ_ID=G.ID.shape[3], NX_T=1, NY_T=1, NZ_T=1), options=compiler_opts)
+        kernel_fields = SourceModule(kernel_template_fields.substitute(REAL=config.dtypes['C_float_or_double'], REAL_OR_COMPLEX=config.dtypes['C_complex'], N_updatecoeffsE=G.updatecoeffsE.size, N_updatecoeffsH=G.updatecoeffsH.size, NY_MATCOEFFS=G.updatecoeffsE.shape[1], NY_MATDISPCOEFFS=1, NX_FIELDS=G.Ex.shape[0], NY_FIELDS=G.Ex.shape[1], NZ_FIELDS=G.Ex.shape[2], NX_ID=G.ID.shape[1], NY_ID=G.ID.shape[2], NZ_ID=G.ID.shape[3], NX_T=1, NY_T=1, NZ_T=1), options=compiler_opts)
 
     else:   # # If there are any dispersive materials (updates are split into two parts as they require present and updated electric field values).
-        kernels_fields = SourceModule(kernels_template_fields.substitute(REAL=config.dtypes['C_float_or_double'], REAL_OR_COMPLEX=config.dtypes['C_complex'], N_updatecoeffsE=G.updatecoeffsE.size, N_updatecoeffsH=G.updatecoeffsH.size, NY_MATCOEFFS=G.updatecoeffsE.shape[1], NY_MATDISPCOEFFS=G.updatecoeffsdispersive.shape[1], NX_FIELDS=G.Ex.shape[0], NY_FIELDS=G.Ex.shape[1], NZ_FIELDS=G.Ex.shape[2], NX_ID=G.ID.shape[1], NY_ID=G.ID.shape[2], NZ_ID=G.ID.shape[3], NX_T=G.Tx.shape[1], NY_T=G.Tx.shape[2], NZ_T=G.Tx.shape[3]), options=compiler_opts)
+        kernel_fields = SourceModule(kernel_template_fields.substitute(REAL=config.dtypes['C_float_or_double'], REAL_OR_COMPLEX=config.dtypes['C_complex'], N_updatecoeffsE=G.updatecoeffsE.size, N_updatecoeffsH=G.updatecoeffsH.size, NY_MATCOEFFS=G.updatecoeffsE.shape[1], NY_MATDISPCOEFFS=G.updatecoeffsdispersive.shape[1], NX_FIELDS=G.Ex.shape[0], NY_FIELDS=G.Ex.shape[1], NZ_FIELDS=G.Ex.shape[2], NX_ID=G.ID.shape[1], NY_ID=G.ID.shape[2], NZ_ID=G.ID.shape[3], NX_T=G.Tx.shape[1], NY_T=G.Tx.shape[2], NZ_T=G.Tx.shape[3]), options=compiler_opts)
-        update_e_dispersive_A_gpu = kernels_fields.get_function("update_e_dispersive_A")
+        update_e_dispersive_A_gpu = kernel_fields.get_function("update_e_dispersive_A")
-        update_e_dispersive_B_gpu = kernels_fields.get_function("update_e_dispersive_B")
+        update_e_dispersive_B_gpu = kernel_fields.get_function("update_e_dispersive_B")
         G.gpu_initialise_dispersive_arrays()
-    update_e_gpu = kernels_fields.get_function("update_e")
+    update_e_gpu = kernel_fields.get_function("update_e")
-    update_h_gpu = kernels_fields.get_function("update_h")
+    update_h_gpu = kernel_fields.get_function("update_h")
 
     # Copy material coefficient arrays to constant memory of GPU (must be <64KB) for fields kernels
-    updatecoeffsE = kernels_fields.get_global('updatecoeffsE')[0]
+    updatecoeffsE = kernel_fields.get_global('updatecoeffsE')[0]
-    updatecoeffsH = kernels_fields.get_global('updatecoeffsH')[0]
+    updatecoeffsH = kernel_fields.get_global('updatecoeffsH')[0]
     if G.updatecoeffsE.nbytes + G.updatecoeffsH.nbytes > config.cuda['gpus'].constmem:
         raise GeneralError('Too many materials in the model to fit onto constant memory of size {} on {} - {} GPU'.format(human_size(config.cuda['gpus'].constmem), config.cuda['gpus'].deviceID, config.cuda['gpus'].name))
     else:
@@ -547,21 +548,21 @@ def solve_gpu(currentmodelrun, modelend, G):
 
     # Sources - initialise arrays on GPU, prepare kernel and get kernel functions
     if G.voltagesources + G.hertziandipoles + G.magneticdipoles:
-        kernels_sources = SourceModule(kernels_template_sources.substitute(REAL=config.dtypes['C_float_or_double'], N_updatecoeffsE=G.updatecoeffsE.size, N_updatecoeffsH=G.updatecoeffsH.size, NY_MATCOEFFS=G.updatecoeffsE.shape[1], NY_SRCINFO=4, NY_SRCWAVES=G.iterations, NX_FIELDS=G.Ex.shape[0], NY_FIELDS=G.Ex.shape[1], NZ_FIELDS=G.Ex.shape[2], NX_ID=G.ID.shape[1], NY_ID=G.ID.shape[2], NZ_ID=G.ID.shape[3]), options=compiler_opts)
+        kernel_sources = SourceModule(kernel_template_sources.substitute(REAL=config.dtypes['C_float_or_double'], N_updatecoeffsE=G.updatecoeffsE.size, N_updatecoeffsH=G.updatecoeffsH.size, NY_MATCOEFFS=G.updatecoeffsE.shape[1], NY_SRCINFO=4, NY_SRCWAVES=G.iterations, NX_FIELDS=G.Ex.shape[0], NY_FIELDS=G.Ex.shape[1], NZ_FIELDS=G.Ex.shape[2], NX_ID=G.ID.shape[1], NY_ID=G.ID.shape[2], NZ_ID=G.ID.shape[3]), options=compiler_opts)
         # Copy material coefficient arrays to constant memory of GPU (must be <64KB) for source kernels
-        updatecoeffsE = kernels_sources.get_global('updatecoeffsE')[0]
+        updatecoeffsE = kernel_sources.get_global('updatecoeffsE')[0]
-        updatecoeffsH = kernels_sources.get_global('updatecoeffsH')[0]
+        updatecoeffsH = kernel_sources.get_global('updatecoeffsH')[0]
         drv.memcpy_htod(updatecoeffsE, G.updatecoeffsE)
         drv.memcpy_htod(updatecoeffsH, G.updatecoeffsH)
         if G.hertziandipoles:
             srcinfo1_hertzian_gpu, srcinfo2_hertzian_gpu, srcwaves_hertzian_gpu = gpu_initialise_src_arrays(G.hertziandipoles, G)
-            update_hertzian_dipole_gpu = kernels_sources.get_function("update_hertzian_dipole")
+            update_hertzian_dipole_gpu = kernel_sources.get_function("update_hertzian_dipole")
         if G.magneticdipoles:
             srcinfo1_magnetic_gpu, srcinfo2_magnetic_gpu, srcwaves_magnetic_gpu = gpu_initialise_src_arrays(G.magneticdipoles, G)
-            update_magnetic_dipole_gpu = kernels_sources.get_function("update_magnetic_dipole")
+            update_magnetic_dipole_gpu = kernel_sources.get_function("update_magnetic_dipole")
         if G.voltagesources:
             srcinfo1_voltage_gpu, srcinfo2_voltage_gpu, srcwaves_voltage_gpu = gpu_initialise_src_arrays(G.voltagesources, G)
-            update_voltage_source_gpu = kernels_sources.get_function("update_voltage_source")
+            update_voltage_source_gpu = kernel_sources.get_function("update_voltage_source")
 
     # Snapshots - initialise arrays on GPU, prepare kernel and get kernel functions
     if G.snapshots:

gprMax/pml.py

@@ -29,11 +29,11 @@ class CFSParameter(object):
 
     # Allowable scaling profiles and directions
     scalingprofiles = {'constant': 0, 'linear': 1, 'quadratic': 2, 'cubic': 3,
-                        'quartic': 4, 'quintic': 5, 'sextic': 6, 'septic': 7, 'octic': 8}
+                       'quartic': 4, 'quintic': 5, 'sextic': 6, 'septic': 7, 'octic': 8}
     scalingdirections = ['forward', 'reverse']
 
     def __init__(self, ID=None, scaling='polynomial', scalingprofile=None,
-                    scalingdirection='forward', min=0, max=0):
+                 scalingdirection='forward', min=0, max=0):
         """
         Args:
             ID (str): Identifier for CFS parameter, can be: 'alpha', 'kappa' or 'sigma'.
@@ -102,7 +102,7 @@ class CFS(object):
         """
 
         tmp = (np.linspace(0, (len(Evalues) - 1) + 0.5, num=2 * len(Evalues))
-                / (len(Evalues) - 1)) ** order
+               / (len(Evalues) - 1)) ** order
         Evalues = tmp[0:-1:2]
         Hvalues = tmp[1::2]
 
@@ -133,8 +133,8 @@ class CFS(object):
 
         elif parameter.scaling == 'polynomial':
             Evalues, Hvalues = self.scaling_polynomial(
-                                CFSParameter.scalingprofiles[parameter.scalingprofile],
+                CFSParameter.scalingprofiles[parameter.scalingprofile],
-                                Evalues, Hvalues)
+                Evalues, Hvalues)
             if parameter.ID == 'alpha':
                 Evalues = Evalues * (self.alpha.max - self.alpha.min) + self.alpha.min
                 Hvalues = Hvalues * (self.alpha.max - self.alpha.min) + self.alpha.min
@@ -163,7 +163,7 @@ class PML(object):
 
     # Available PML formulations:
     # Higher Order RIPML (HORIPML) see: https://doi.org/10.1109/TAP.2011.2180344
-    # Multipole RIPML (MRIPML) see: https://doi.org/10.1109/TAP.2018.2823864
+    # Multipole RIPML (MRIPML) see: https://doi.org/10.1109/TAP.2018.2823864
     formulations = ['HORIPML', 'MRIPML']
 
     # PML slabs IDs at boundaries of domain.
@@ -216,31 +216,31 @@ class PML(object):
 
         if self.direction[0] == 'x':
             self.EPhi1 = np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz + 1),
-                            dtype=config.dtypes['float_or_double'])
+                                  dtype=config.dtypes['float_or_double'])
             self.EPhi2 = np.zeros((len(self.CFS), self.nx + 1, self.ny + 1, self.nz),
-                            dtype=config.dtypes['float_or_double'])
+                                  dtype=config.dtypes['float_or_double'])
             self.HPhi1 = np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz),
-                            dtype=config.dtypes['float_or_double'])
+                                  dtype=config.dtypes['float_or_double'])
             self.HPhi2 = np.zeros((len(self.CFS), self.nx, self.ny, self.nz + 1),
-                            dtype=config.dtypes['float_or_double'])
+                                  dtype=config.dtypes['float_or_double'])
         elif self.direction[0] == 'y':
             self.EPhi1 = np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz + 1),
-                            dtype=config.dtypes['float_or_double'])
+                                  dtype=config.dtypes['float_or_double'])
             self.EPhi2 = np.zeros((len(self.CFS), self.nx + 1, self.ny + 1, self.nz),
-                            dtype=config.dtypes['float_or_double'])
+                                  dtype=config.dtypes['float_or_double'])
             self.HPhi1 = np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz),
-                            dtype=config.dtypes['float_or_double'])
+                                  dtype=config.dtypes['float_or_double'])
             self.HPhi2 = np.zeros((len(self.CFS), self.nx, self.ny, self.nz + 1),
-                            dtype=config.dtypes['float_or_double'])
+                                  dtype=config.dtypes['float_or_double'])
         elif self.direction[0] == 'z':
             self.EPhi1 = np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz + 1),
-                            dtype=config.dtypes['float_or_double'])
+                                  dtype=config.dtypes['float_or_double'])
             self.EPhi2 = np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz + 1),
-                            dtype=config.dtypes['float_or_double'])
+                                  dtype=config.dtypes['float_or_double'])
             self.HPhi1 = np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz),
-                            dtype=config.dtypes['float_or_double'])
+                                  dtype=config.dtypes['float_or_double'])
             self.HPhi2 = np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz),
-                            dtype=config.dtypes['float_or_double'])
+                                  dtype=config.dtypes['float_or_double'])
 
     def calculate_update_coeffs(self, er, mr, G):
         """Calculates electric and magnetic update coefficients for the PML.
@@ -253,21 +253,21 @@ class PML(object):
         """
 
         self.ERA = np.zeros((len(self.CFS), self.thickness),
-                    dtype=config.dtypes['float_or_double'])
+                            dtype=config.dtypes['float_or_double'])
         self.ERB = np.zeros((len(self.CFS), self.thickness),
-                    dtype=config.dtypes['float_or_double'])
+                            dtype=config.dtypes['float_or_double'])
         self.ERE = np.zeros((len(self.CFS), self.thickness),
-                    dtype=config.dtypes['float_or_double'])
+                            dtype=config.dtypes['float_or_double'])
         self.ERF = np.zeros((len(self.CFS), self.thickness),
-                    dtype=config.dtypes['float_or_double'])
+                            dtype=config.dtypes['float_or_double'])
         self.HRA = np.zeros((len(self.CFS), self.thickness),
-                    dtype=config.dtypes['float_or_double'])
+                            dtype=config.dtypes['float_or_double'])
         self.HRB = np.zeros((len(self.CFS), self.thickness),
-                    dtype=config.dtypes['float_or_double'])
+                            dtype=config.dtypes['float_or_double'])
         self.HRE = np.zeros((len(self.CFS), self.thickness),
-                    dtype=config.dtypes['float_or_double'])
+                            dtype=config.dtypes['float_or_double'])
         self.HRF = np.zeros((len(self.CFS), self.thickness),
-                    dtype=config.dtypes['float_or_double'])
+                            dtype=config.dtypes['float_or_double'])
 
         for x, cfs in enumerate(self.CFS):
             if not cfs.sigma.max:
@@ -283,7 +283,7 @@ class PML(object):
                 self.ERA[x, :] = (2 * config.e0 + G.dt * Ealpha) / tmp
                 self.ERB[x, :] = (2 * config.e0 * Ekappa) / tmp
                 self.ERE[x, :] = ((2 * config.e0 * Ekappa) - G.dt
-                                    * (Ealpha * Ekappa + Esigma)) / tmp
+                                  * (Ealpha * Ekappa + Esigma)) / tmp
                 self.ERF[x, :] = (2 * Esigma * G.dt) / (Ekappa * tmp)
 
                 # HORIPML magnetic update coefficients
@@ -291,7 +291,7 @@ class PML(object):
                 self.HRA[x, :] = (2 * config.e0 + G.dt * Halpha) / tmp
                 self.HRB[x, :] = (2 * config.e0 * Hkappa) / tmp
                 self.HRE[x, :] = ((2 * config.e0 * Hkappa) - G.dt
-                                    * (Halpha * Hkappa + Hsigma)) / tmp
+                                  * (Halpha * Hkappa + Hsigma)) / tmp
                 self.HRF[x, :] = (2 * Hsigma * G.dt) / (Hkappa * tmp)
 
             elif G.pmlformulation == 'MRIPML':
@@ -318,9 +318,9 @@ class PML(object):
         pmlmodule = 'gprMax.cython.pml_updates_electric_' + G.pmlformulation
         func = getattr(import_module(pmlmodule), 'order' + str(len(self.CFS)) + '_' + self.direction)
         func(self.xs, self.xf, self.ys, self.yf, self.zs, self.zf,
-                config.hostinfo['ompthreads'], G.updatecoeffsE, G.ID,
+             config.hostinfo['ompthreads'], G.updatecoeffsE, G.ID,
-                G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, self.EPhi1, self.EPhi2,
+             G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, self.EPhi1, self.EPhi2,
-                self.ERA, self.ERB, self.ERE, self.ERF, self.d)
+             self.ERA, self.ERB, self.ERE, self.ERF, self.d)
 
     def update_magnetic(self, G):
         """This functions updates magnetic field components with the PML correction.
@@ -332,9 +332,9 @@ class PML(object):
         pmlmodule = 'gprMax.cython.pml_updates_magnetic_' + G.pmlformulation
         func = getattr(import_module(pmlmodule), 'order' + str(len(self.CFS)) + '_' + self.direction)
         func(self.xs, self.xf, self.ys, self.yf, self.zs, self.zf,
-                config.hostinfo['ompthreads'], G.updatecoeffsH, G.ID,
+             config.hostinfo['ompthreads'], G.updatecoeffsH, G.ID,
-                G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, self.HPhi1, self.HPhi2,
+             G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, self.HPhi1, self.HPhi2,
-                self.HRA, self.HRB, self.HRE, self.HRF, self.d)
+             self.HRA, self.HRB, self.HRE, self.HRF, self.d)
 
     def gpu_set_blocks_per_grid(self, G):
         """Set the blocks per grid size used for updating the PML field arrays on a GPU.
@@ -344,9 +344,9 @@ class PML(object):
         """
 
         config.cuda['gpus'].bpg = (int(np.ceil(((self.EPhi1.shape[1] + 1)
-                                    * (self.EPhi1.shape[2] + 1)
+                                                * (self.EPhi1.shape[2] + 1)
-                                    * (self.EPhi1.shape[3] + 1))
+                                                * (self.EPhi1.shape[3] + 1))
-                                    / config.cuda['gpus'].tpb[0])), 1, 1)
+                                               / config.cuda['gpus'].tpb[0])), 1, 1)
 
     def gpu_initialise_arrays(self):
         """Initialise PML field and coefficient arrays on GPU."""
@@ -379,9 +379,9 @@ class PML(object):
         from pycuda.compiler import SourceModule
 
         self.update_electric = kernelselectric.get_function('order' + str(len(self.CFS))
-                                + '_' + self.direction)
+                                                            + '_' + self.direction)
         self.update_magnetic = kernelsmagnetic.get_function('order' + str(len(self.CFS))
-                                + '_' + self.direction)
+                                                            + '_' + self.direction)
 
     def gpu_update_electric(self, G):
         """This functions updates electric field components with the PML
@@ -393,19 +393,19 @@ class PML(object):
         """
 
         self.update_electric(np.int32(self.xs), np.int32(self.xf),
-                                np.int32(self.ys), np.int32(self.yf),
+                             np.int32(self.ys), np.int32(self.yf),
-                                np.int32(self.zs), np.int32(self.zf),
+                             np.int32(self.zs), np.int32(self.zf),
-                                np.int32(self.EPhi1.shape[1]), np.int32(self.EPhi1.shape[2]),
+                             np.int32(self.EPhi1.shape[1]), np.int32(self.EPhi1.shape[2]),
-                                np.int32(self.EPhi1.shape[3]), np.int32(self.EPhi2.shape[1]),
+                             np.int32(self.EPhi1.shape[3]), np.int32(self.EPhi2.shape[1]),
-                                np.int32(self.EPhi2.shape[2]), np.int32(self.EPhi2.shape[3]),
+                             np.int32(self.EPhi2.shape[2]), np.int32(self.EPhi2.shape[3]),
-                                np.int32(self.thickness), G.ID_gpu.gpudata,
+                             np.int32(self.thickness), G.ID_gpu.gpudata,
-                                G.Ex_gpu.gpudata, G.Ey_gpu.gpudata, G.Ez_gpu.gpudata,
+                             G.Ex_gpu.gpudata, G.Ey_gpu.gpudata, G.Ez_gpu.gpudata,
-                                G.Hx_gpu.gpudata, G.Hy_gpu.gpudata, G.Hz_gpu.gpudata,
+                             G.Hx_gpu.gpudata, G.Hy_gpu.gpudata, G.Hz_gpu.gpudata,
-                                self.EPhi1_gpu.gpudata, self.EPhi2_gpu.gpudata,
+                             self.EPhi1_gpu.gpudata, self.EPhi2_gpu.gpudata,
-                                self.ERA_gpu.gpudata, self.ERB_gpu.gpudata,
+                             self.ERA_gpu.gpudata, self.ERB_gpu.gpudata,
-                                self.ERE_gpu.gpudata, self.ERF_gpu.gpudata,
+                             self.ERE_gpu.gpudata, self.ERF_gpu.gpudata,
-                                config.dtypes['float_or_double'](self.d),
+                             config.dtypes['float_or_double'](self.d),
-                                block=config.cuda['gpus'].tpb, grid=config.cuda['gpus'].bpg)
+                             block=config.cuda['gpus'].tpb, grid=config.cuda['gpus'].bpg)
 
     def gpu_update_magnetic(self, G):
         """This functions updates magnetic field components with the PML
@@ -417,19 +417,19 @@ class PML(object):
         """
 
         self.update_magnetic(np.int32(self.xs), np.int32(self.xf),
-                                np.int32(self.ys), np.int32(self.yf),
+                             np.int32(self.ys), np.int32(self.yf),
-                                np.int32(self.zs), np.int32(self.zf),
+                             np.int32(self.zs), np.int32(self.zf),
-                                np.int32(self.HPhi1.shape[1]), np.int32(self.HPhi1.shape[2]),
+                             np.int32(self.HPhi1.shape[1]), np.int32(self.HPhi1.shape[2]),
-                                np.int32(self.HPhi1.shape[3]), np.int32(self.HPhi2.shape[1]),
+                             np.int32(self.HPhi1.shape[3]), np.int32(self.HPhi2.shape[1]),
-                                np.int32(self.HPhi2.shape[2]), np.int32(self.HPhi2.shape[3]),
+                             np.int32(self.HPhi2.shape[2]), np.int32(self.HPhi2.shape[3]),
-                                np.int32(self.thickness), G.ID_gpu.gpudata,
+                             np.int32(self.thickness), G.ID_gpu.gpudata,
-                                G.Ex_gpu.gpudata, G.Ey_gpu.gpudata, G.Ez_gpu.gpudata,
+                             G.Ex_gpu.gpudata, G.Ey_gpu.gpudata, G.Ez_gpu.gpudata,
-                                G.Hx_gpu.gpudata, G.Hy_gpu.gpudata, G.Hz_gpu.gpudata,
+                             G.Hx_gpu.gpudata, G.Hy_gpu.gpudata, G.Hz_gpu.gpudata,
-                                self.HPhi1_gpu.gpudata, self.HPhi2_gpu.gpudata,
+                             self.HPhi1_gpu.gpudata, self.HPhi2_gpu.gpudata,
-                                self.HRA_gpu.gpudata, self.HRB_gpu.gpudata,
+                             self.HRA_gpu.gpudata, self.HRB_gpu.gpudata,
-                                self.HRE_gpu.gpudata, self.HRF_gpu.gpudata,
+                             self.HRE_gpu.gpudata, self.HRF_gpu.gpudata,
-                                config.dtypes['float_or_double'](self.d),
+                             config.dtypes['float_or_double'](self.d),
-                                block=config.cuda['gpus'].tpb, grid=config.cuda['gpus'].bpg)
+                             block=config.cuda['gpus'].tpb, grid=config.cuda['gpus'].bpg)
 
 
 def build_pmls(G, pbar):

gprMax/receivers.py

@@ -64,7 +64,7 @@ def gpu_initialise_rx_arrays(G):
     # Array to store field components for receivers on GPU - rows are field components;
     # columns are iterations; pages are receivers
     rxs = np.zeros((Rx.maxnumoutputs, G.iterations, len(G.rxs)),
-            dtype=config.dtypes['float_or_double'])
+                   dtype=config.dtypes['float_or_double'])
 
     # Copy arrays to GPU
     rxcoords_gpu = gpuarray.to_gpu(rxcoords)
@@ -89,5 +89,5 @@ def gpu_get_rx_array(rxs_gpu, rxcoords_gpu, G):
             if rx.xcoord == rxcoords_gpu[rxgpu, 0] and \
                rx.ycoord == rxcoords_gpu[rxgpu, 1] and \
                rx.zcoord == rxcoords_gpu[rxgpu, 2]:
-                for k, v in rx.outputs.items():
+                for k in rx.outputs.items():
-                    v = rxs_gpu[Rx.gpu_allowableoutputs.index(k), :, rxgpu]
+                    rx.outputs[k] = rxs_gpu[Rx.gpu_allowableoutputs.index(k), :, rxgpu]

gprMax/snapshots.py

@@ -55,7 +55,7 @@ class Snapshot(object):
         floatstring = 'd'
 
     def __init__(self, xs=None, ys=None, zs=None, xf=None, yf=None, zf=None,
-                dx=None, dy=None, dz=None, time=None, filename=None):
+                 dx=None, dy=None, dz=None, time=None, filename=None):
         """
         Args:
             xs, xf, ys, yf, zs, zf (int): Extent of the volume in cells.
@@ -84,10 +84,10 @@ class Snapshot(object):
         self.datasizefield = (3 * np.dtype(config.dtypes['float_or_double']).itemsize
                                 * self.ncells)
         self.vtkdatawritesize = ((self.fieldoutputs['electric']
-                                + self.fieldoutputs['magnetic']) * self.datasizefield
+                                  + self.fieldoutputs['magnetic']) * self.datasizefield
-                                + (self.fieldoutputs['electric']
+                                 + (self.fieldoutputs['electric']
-                                + self.fieldoutputs['magnetic'])
+                                    + self.fieldoutputs['magnetic'])
-                                * np.dtype(np.uint32).itemsize)
+                                 * np.dtype(np.uint32).itemsize)
         self.time = time
         self.basefilename = filename
 
@@ -147,36 +147,36 @@ class Snapshot(object):
         """
 
         hfield_offset = (3 * np.dtype(config.dtypes['float_or_double']).itemsize
-                        * self.ncells + np.dtype(np.uint32).itemsize)
+                         * self.ncells + np.dtype(np.uint32).itemsize)
 
         self.filehandle = open(self.filename, 'wb')
         self.filehandle.write('<?xml version="1.0"?>\n'.encode('utf-8'))
         self.filehandle.write('<VTKFile type="ImageData" version="1.0" byte_order="{}">\n'
-                                .format(Snapshot.byteorder).encode('utf-8'))
+                              .format(Snapshot.byteorder).encode('utf-8'))
         self.filehandle.write('<ImageData WholeExtent="{} {} {} {} {} {}" Origin="0 0 0" Spacing="{:.3} {:.3} {:.3}">\n'
-                                .format(self.xs, round_value(self.xf / self.dx),
+                              .format(self.xs, round_value(self.xf / self.dx),
-                                self.ys, round_value(self.yf / self.dy), self.zs,
+                                      self.ys, round_value(self.yf / self.dy), self.zs,
-                                round_value(self.zf / self.dz), self.dx * G.dx,
+                                      round_value(self.zf / self.dz), self.dx * G.dx,
-                                self.dy * G.dy, self.dz * G.dz).encode('utf-8'))
+                                      self.dy * G.dy, self.dz * G.dz).encode('utf-8'))
         self.filehandle.write('<Piece Extent="{} {} {} {} {} {}">\n'
-                                .format(self.xs, round_value(self.xf / self.dx),
+                              .format(self.xs, round_value(self.xf / self.dx),
-                                self.ys, round_value(self.yf / self.dy),
+                                      self.ys, round_value(self.yf / self.dy),
-                                self.zs, round_value(self.zf / self.dz)).encode('utf-8'))
+                                      self.zs, round_value(self.zf / self.dz)).encode('utf-8'))
 
         if self.fieldoutputs['electric'] and self.fieldoutputs['magnetic']:
             self.filehandle.write('<CellData Vectors="E-field H-field">\n'.encode('utf-8'))
             self.filehandle.write('<DataArray type="{}" Name="E-field" NumberOfComponents="3" format="appended" offset="0" />\n'
-                                    .format(Snapshot.floatname).encode('utf-8'))
+                                  .format(Snapshot.floatname).encode('utf-8'))
             self.filehandle.write('<DataArray type="{}" Name="H-field" NumberOfComponents="3" format="appended" offset="{}" />\n'
-                                    .format(Snapshot.floatname, hfield_offset).encode('utf-8'))
+                                  .format(Snapshot.floatname, hfield_offset).encode('utf-8'))
         elif self.fieldoutputs['electric']:
             self.filehandle.write('<CellData Vectors="E-field">\n'.encode('utf-8'))
             self.filehandle.write('<DataArray type="{}" Name="E-field" NumberOfComponents="3" format="appended" offset="0" />\n'
-                                    .format(Snapshot.floatname).encode('utf-8'))
+                                  .format(Snapshot.floatname).encode('utf-8'))
         elif self.fieldoutputs['magnetic']:
             self.filehandle.write('<CellData Vectors="H-field">\n'.encode('utf-8'))
             self.filehandle.write('<DataArray type="{}" Name="H-field" NumberOfComponents="3" format="appended" offset="0" />\n'
-                                    .format(Snapshot.floatname).encode('utf-8'))
+                                  .format(Snapshot.floatname).encode('utf-8'))
 
         self.filehandle.write('</CellData>\n</Piece>\n</ImageData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
 
@@ -224,20 +224,20 @@ def gpu_initialise_snapshot_array(G):
 
     # 4D arrays to store snapshots on GPU, e.g. snapEx(time, x, y, z);
     # if snapshots are not being stored on the GPU during the simulation then
-    # they are copied back to the host after each iteration, hence numsnaps = 1
+    # they are copied back to the host after each iteration, hence numsnaps = 1
     numsnaps = 1 if config.cuda['snapsgpu2cpu'] else len(G.snapshots)
     snapEx = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
-                        dtype=config.dtypes['float_or_double'])
+                      dtype=config.dtypes['float_or_double'])
     snapEy = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
-                        dtype=config.dtypes['float_or_double'])
+                      dtype=config.dtypes['float_or_double'])
     snapEz = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
-                        dtype=config.dtypes['float_or_double'])
+                      dtype=config.dtypes['float_or_double'])
     snapHx = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
-                        dtype=config.dtypes['float_or_double'])
+                      dtype=config.dtypes['float_or_double'])
     snapHy = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
-                        dtype=config.dtypes['float_or_double'])
+                      dtype=config.dtypes['float_or_double'])
     snapHz = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
-                        dtype=config.dtypes['float_or_double'])
+                      dtype=config.dtypes['float_or_double'])
 
     # Copy arrays to GPU
     snapEx_gpu = gpuarray.to_gpu(snapEx)
@@ -260,8 +260,8 @@ def gpu_get_snapshot_array(snapEx_gpu, snapEy_gpu, snapEz_gpu, snapHx_gpu, snapH
     """
 
     snap.electric = np.stack((snapEx_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf],
-                                snapEy_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf],
+                              snapEy_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf],
-                                snapEz_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf])).reshape(-1, order='F')
+                              snapEz_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf])).reshape(-1, order='F')
     snap.magnetic = np.stack((snapHx_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf],
-                                snapHy_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf],
+                              snapHy_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf],
-                                snapHz_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf])).reshape(-1, order='F')
+                              snapHz_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf])).reshape(-1, order='F')

gprMax/utilities.py

@@ -414,6 +414,7 @@ def detect_check_gpus(deviceIDs):
 
     return gpus, allgpustext
 
+
 def timer():
     """Function to return the current process wide time in fractional seconds."""
     return perf_counter()