diff --git a/gprMax/config.py b/gprMax/config.py index 2809c3d0..1e310152 100644 --- a/gprMax/config.py +++ b/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, - 'C_float_or_double': 'float', 'C_complex': 'pycuda::complex'} + 'cython_float_or_double': cython.float, 'cython_complex': cython.floatcomplex, + 'C_float_or_double': 'float', 'C_complex': 'pycuda::complex'} elif precision == 'double': dtypes = {'float_or_double': np.float64, 'complex': np.complex128, - 'cython_float_or_double': cython.double, 'cython_complex': cython.doublecomplex, - 'C_float_or_double': 'double', 'C_complex': 'pycuda::complex'} + 'cython_float_or_double': cython.double, 'cython_complex': cython.doublecomplex, + 'C_float_or_double': 'double', 'C_complex': 'pycuda::complex'} diff --git a/gprMax/fields_outputs.py b/gprMax/fields_outputs.py index 2a89acd9..57b8f0b1 100644 --- a/gprMax/fields_outputs.py +++ b/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] diff --git a/gprMax/fractals.py b/gprMax/fractals.py index 07d8f7e1..3822e016 100644 --- a/gprMax/fractals.py +++ b/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] - ((self.fractalrange[1] - - self.fractalrange[0]) / fractalrange) * fractalmin) + - self.fractalrange[0]) / fractalrange) + + self.fractalrange[0] - ((self.fractalrange[1] + - 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) diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index c02e898b..b71d1d94 100644 --- a/gprMax/gprMax.py +++ b/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) diff --git a/gprMax/input_cmds_file.py b/gprMax/input_cmds_file.py index 4c39e73d..b2d41f69 100644 --- a/gprMax/input_cmds_file.py +++ b/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: diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py index 2a6ad7ae..da97d82c 100644 --- a/gprMax/input_cmds_multiuse.py +++ b/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 diff --git a/gprMax/materials.py b/gprMax/materials.py index 0a82a5bc..c9097f61 100644 --- a/gprMax/materials.py +++ b/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]', - 'gamma\n[Hz]', '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']] + 'Delta\neps_r', 'tau\n[s]', 'omega\n[Hz]', 'delta\n[Hz]', + '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 diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py index 09375e26..bd6683a7 100644 --- a/gprMax/model_build_run.py +++ b/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) - update_e_dispersive_A_gpu = kernels_fields.get_function("update_e_dispersive_A") - update_e_dispersive_B_gpu = kernels_fields.get_function("update_e_dispersive_B") + 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 = kernel_fields.get_function("update_e_dispersive_A") + 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_h_gpu = kernels_fields.get_function("update_h") + update_e_gpu = kernel_fields.get_function("update_e") + 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] - updatecoeffsH = kernels_fields.get_global('updatecoeffsH')[0] + updatecoeffsE = kernel_fields.get_global('updatecoeffsE')[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] - updatecoeffsH = kernels_sources.get_global('updatecoeffsH')[0] + updatecoeffsE = kernel_sources.get_global('updatecoeffsE')[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: diff --git a/gprMax/pml.py b/gprMax/pml.py index 488b67a6..59547254 100644 --- a/gprMax/pml.py +++ b/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], - Evalues, Hvalues) + CFSParameter.scalingprofiles[parameter.scalingprofile], + 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, - 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) + config.hostinfo['ompthreads'], G.updatecoeffsE, G.ID, + 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) 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, - 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) + config.hostinfo['ompthreads'], G.updatecoeffsH, G.ID, + 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) 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[3] + 1)) - / config.cuda['gpus'].tpb[0])), 1, 1) + * (self.EPhi1.shape[2] + 1) + * (self.EPhi1.shape[3] + 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.zs), np.int32(self.zf), - 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.EPhi2.shape[2]), np.int32(self.EPhi2.shape[3]), - np.int32(self.thickness), G.ID_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, - self.EPhi1_gpu.gpudata, self.EPhi2_gpu.gpudata, - self.ERA_gpu.gpudata, self.ERB_gpu.gpudata, - self.ERE_gpu.gpudata, self.ERF_gpu.gpudata, - config.dtypes['float_or_double'](self.d), - block=config.cuda['gpus'].tpb, grid=config.cuda['gpus'].bpg) + np.int32(self.ys), np.int32(self.yf), + 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[3]), np.int32(self.EPhi2.shape[1]), + np.int32(self.EPhi2.shape[2]), np.int32(self.EPhi2.shape[3]), + np.int32(self.thickness), G.ID_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, + self.EPhi1_gpu.gpudata, self.EPhi2_gpu.gpudata, + self.ERA_gpu.gpudata, self.ERB_gpu.gpudata, + self.ERE_gpu.gpudata, self.ERF_gpu.gpudata, + config.dtypes['float_or_double'](self.d), + 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.zs), np.int32(self.zf), - 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.HPhi2.shape[2]), np.int32(self.HPhi2.shape[3]), - np.int32(self.thickness), G.ID_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, - self.HPhi1_gpu.gpudata, self.HPhi2_gpu.gpudata, - self.HRA_gpu.gpudata, self.HRB_gpu.gpudata, - self.HRE_gpu.gpudata, self.HRF_gpu.gpudata, - config.dtypes['float_or_double'](self.d), - block=config.cuda['gpus'].tpb, grid=config.cuda['gpus'].bpg) + np.int32(self.ys), np.int32(self.yf), + 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[3]), np.int32(self.HPhi2.shape[1]), + np.int32(self.HPhi2.shape[2]), np.int32(self.HPhi2.shape[3]), + np.int32(self.thickness), G.ID_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, + self.HPhi1_gpu.gpudata, self.HPhi2_gpu.gpudata, + self.HRA_gpu.gpudata, self.HRB_gpu.gpudata, + self.HRE_gpu.gpudata, self.HRF_gpu.gpudata, + config.dtypes['float_or_double'](self.d), + block=config.cuda['gpus'].tpb, grid=config.cuda['gpus'].bpg) def build_pmls(G, pbar): diff --git a/gprMax/receivers.py b/gprMax/receivers.py index 861df298..e7e10bb4 100644 --- a/gprMax/receivers.py +++ b/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(): - v = rxs_gpu[Rx.gpu_allowableoutputs.index(k), :, rxgpu] + for k in rx.outputs.items(): + rx.outputs[k] = rxs_gpu[Rx.gpu_allowableoutputs.index(k), :, rxgpu] diff --git a/gprMax/snapshots.py b/gprMax/snapshots.py index 9b02056f..1caf678d 100644 --- a/gprMax/snapshots.py +++ b/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['electric'] - + self.fieldoutputs['magnetic']) - * np.dtype(np.uint32).itemsize) + + self.fieldoutputs['magnetic']) * self.datasizefield + + (self.fieldoutputs['electric'] + + self.fieldoutputs['magnetic']) + * 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('\n'.encode('utf-8')) self.filehandle.write('\n' - .format(Snapshot.byteorder).encode('utf-8')) + .format(Snapshot.byteorder).encode('utf-8')) self.filehandle.write('\n' - .format(self.xs, round_value(self.xf / self.dx), - self.ys, round_value(self.yf / self.dy), self.zs, - round_value(self.zf / self.dz), self.dx * G.dx, - self.dy * G.dy, self.dz * G.dz).encode('utf-8')) + .format(self.xs, round_value(self.xf / self.dx), + self.ys, round_value(self.yf / self.dy), self.zs, + round_value(self.zf / self.dz), self.dx * G.dx, + self.dy * G.dy, self.dz * G.dz).encode('utf-8')) self.filehandle.write('\n' - .format(self.xs, round_value(self.xf / self.dx), - self.ys, round_value(self.yf / self.dy), - self.zs, round_value(self.zf / self.dz)).encode('utf-8')) + .format(self.xs, round_value(self.xf / self.dx), + self.ys, round_value(self.yf / self.dy), + self.zs, round_value(self.zf / self.dz)).encode('utf-8')) if self.fieldoutputs['electric'] and self.fieldoutputs['magnetic']: self.filehandle.write('\n'.encode('utf-8')) self.filehandle.write('\n' - .format(Snapshot.floatname).encode('utf-8')) + .format(Snapshot.floatname).encode('utf-8')) self.filehandle.write('\n' - .format(Snapshot.floatname, hfield_offset).encode('utf-8')) + .format(Snapshot.floatname, hfield_offset).encode('utf-8')) elif self.fieldoutputs['electric']: self.filehandle.write('\n'.encode('utf-8')) self.filehandle.write('\n' - .format(Snapshot.floatname).encode('utf-8')) + .format(Snapshot.floatname).encode('utf-8')) elif self.fieldoutputs['magnetic']: self.filehandle.write('\n'.encode('utf-8')) self.filehandle.write('\n' - .format(Snapshot.floatname).encode('utf-8')) + .format(Snapshot.floatname).encode('utf-8')) self.filehandle.write('\n\n\n\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], - snapEz_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf])).reshape(-1, order='F') + 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') 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], - snapHz_gpu[i, snap.xs:snap.xf, snap.ys:snap.yf, snap.zs:snap.zf])).reshape(-1, order='F') + 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') diff --git a/gprMax/utilities.py b/gprMax/utilities.py index 9ebd682d..d9a8cc30 100644 --- a/gprMax/utilities.py +++ b/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()