diff --git a/gprMax/__init__.py b/gprMax/__init__.py index 92526be3..15cfe14e 100644 --- a/gprMax/__init__.py +++ b/gprMax/__init__.py @@ -50,7 +50,7 @@ from .cmds_geometry.add_surface_roughness import AddSurfaceRoughness from .cmds_geometry.add_surface_water import AddSurfaceWater from .cmds_geometry.add_grass import AddGrass -from .input_cmds_file import user_libs_fn_to_scene_obj +from .hash_cmds_file import user_libs_fn_to_scene_obj from .scene import Scene from .gprMax import run as run diff --git a/gprMax/cmds_multiple.py b/gprMax/cmds_multiple.py index 7fe1624a..adecadf4 100644 --- a/gprMax/cmds_multiple.py +++ b/gprMax/cmds_multiple.py @@ -39,9 +39,9 @@ from .subgrids.base import SubGridBase from .utilities import round_value from .waveforms import Waveform as WaveformUser - log = logging.getLogger(__name__) + class UserObjectMulti: """Object that can occur multiple times in a model.""" @@ -143,26 +143,26 @@ class VoltageSource(UserObjectMulti): resistance = self.kwargs['resistance'] waveform_id = self.kwargs['waveform_id'] except KeyError: - raise CmdInputError(f"'{self.__str__()}' requires at least six parameters") + raise CmdInputError(f"'{self.params_str()}' requires at least six parameters") # Check polarity & position parameters if polarisation not in ('x', 'y', 'z'): - raise CmdInputError(f"'{self.__str__()}' polarisation must be x, y, or z") + raise CmdInputError(f"'{self.params_str()}' polarisation must be x, y, or z") if '2D TMx' in grid.mode and (polarisation == 'y' or polarisation == 'z'): - raise CmdInputError(f"'{self.__str__()}' polarisation must be x in 2D TMx mode") + raise CmdInputError(f"'{self.params_str()}' polarisation must be x in 2D TMx mode") elif '2D TMy' in grid.mode and (polarisation == 'x' or polarisation == 'z'): - raise CmdInputError(f"'{self.__str__()}' polarisation must be y in 2D TMy mode") + raise CmdInputError(f"'{self.params_str()}' polarisation must be y in 2D TMy mode") elif '2D TMz' in grid.mode and (polarisation == 'x' or polarisation == 'y'): - raise CmdInputError(f"'{self.__str__()}' polarisation must be z in 2D TMz mode") + raise CmdInputError(f"'{self.params_str()}' polarisation must be z in 2D TMz mode") - xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.__str__()) + xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.params_str()) if resistance < 0: - raise CmdInputError(f"'{self.__str__()}' requires a source resistance of zero or greater") + raise CmdInputError(f"'{self.params_str()}' requires a source resistance of zero or greater") # Check if there is a waveformID in the waveforms list if not any(x.ID == waveform_id for x in grid.waveforms): - raise CmdInputError(f"'{self.__str__()}' there is no waveform with the identifier {tmp[5]}") + raise CmdInputError(f"'{self.params_str()}' there is no waveform with the identifier {tmp[5]}") v = VoltageSourceUser() v.polarisation = polarisation @@ -178,11 +178,11 @@ class VoltageSource(UserObjectMulti): stop = self.kwargs['stop'] # Check source start & source remove time parameters if start < 0: - raise CmdInputError(f"'{self.__str__()}' delay of the initiation of the source should not be less than zero") + raise CmdInputError(f"'{self.params_str()}' delay of the initiation of the source should not be less than zero") if stop < 0: - raise CmdInputError(f"'{self.__str__()}' time to remove the source should not be less than zero") + raise CmdInputError(f"'{self.params_str()}' time to remove the source should not be less than zero") if stop - start <= 0: - raise CmdInputError(f"'{self.__str__()}' duration of the source should not be zero or less") + raise CmdInputError(f"'{self.params_str()}' duration of the source should not be zero or less") v.start = start if stop > grid.timewindow: v.stop = grid.timewindow @@ -228,23 +228,23 @@ class HertzianDipole(UserObjectMulti): p1 = self.kwargs['p1'] waveform_id = self.kwargs['waveform_id'] except KeyError: - raise CmdInputError(f"'{self.__str__()}' requires at least 3 parameters") + raise CmdInputError(f"'{self.params_str()}' requires at least 3 parameters") # Check polarity & position parameters if polarisation not in ('x', 'y', 'z'): - raise CmdInputError(f"'{self.__str__()}' polarisation must be x, y, or z") + raise CmdInputError(f"'{self.params_str()}' polarisation must be x, y, or z") if '2D TMx' in grid.mode and (polarisation == 'y' or polarisation == 'z'): - raise CmdInputError(f"'{self.__str__()}' polarisation must be x in 2D TMx mode") + raise CmdInputError(f"'{self.params_str()}' polarisation must be x in 2D TMx mode") elif '2D TMy' in grid.mode and (polarisation == 'x' or polarisation == 'z'): - raise CmdInputError(f"'{self.__str__()}' polarisation must be y in 2D TMy mode") + raise CmdInputError(f"'{self.params_str()}' polarisation must be y in 2D TMy mode") elif '2D TMz' in grid.mode and (polarisation == 'x' or polarisation == 'y'): - raise CmdInputError(f"'{self.__str__()}' polarisation must be z in 2D TMz mode") + raise CmdInputError(f"'{self.params_str()}' polarisation must be z in 2D TMz mode") - xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.__str__()) + xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.params_str()) # Check if there is a waveformID in the waveforms list if not any(x.ID == waveform_id for x in grid.waveforms): - raise CmdInputError(f"'{self.__str__()}' there is no waveform with the identifier {tmp[4]}') + raise CmdInputError(f"'{self.params_str()}' there is no waveform with the identifier {tmp[4]}") h = HertzianDipoleUser() h.polarisation = polarisation @@ -271,11 +271,11 @@ class HertzianDipole(UserObjectMulti): start = self.kwargs['start'] stop = self.kwargs['stop'] if start < 0: - raise CmdInputError(f"'{self.__str__()}' delay of the initiation of the source should not be less than zero") + raise CmdInputError(f"'{self.params_str()}' delay of the initiation of the source should not be less than zero") if stop < 0: - raise CmdInputError(f"'{self.__str__()self.__str__()}' time to remove the source should not be less than zero") + raise CmdInputError(f"'{self.params_str()}' time to remove the source should not be less than zero") if stop - start <= 0: - raise CmdInputError(f"'{self.__str__()}' duration of the source should not be zero or less") + raise CmdInputError(f"'{self.params_str()}' duration of the source should not be zero or less") h.start = start if stop > grid.timewindow: h.stop = grid.timewindow @@ -324,23 +324,23 @@ class MagneticDipole(UserObjectMulti): p1 = self.kwargs['p1'] waveform_id = self.kwargs['waveform_id'] except KeyError: - raise CmdInputError(f"'{self.__str__()}' requires at least five parameters") + raise CmdInputError(f"'{self.params_str()}' requires at least five parameters") # Check polarity & position parameters if polarisation not in ('x', 'y', 'z'): - raise CmdInputError(f"'{self.__str__()}' polarisation must be x, y, or z") + raise CmdInputError(f"'{self.params_str()}' polarisation must be x, y, or z") if '2D TMx' in grid.mode and (polarisation == 'y' or polarisation == 'z'): - raise CmdInputError(f"'{self.__str__()}' polarisation must be x in 2D TMx mode") + raise CmdInputError(f"'{self.params_str()}' polarisation must be x in 2D TMx mode") elif '2D TMy' in grid.mode and (polarisation == 'x' or polarisation == 'z'): - raise CmdInputError(f"'{self.__str__()}' polarisation must be y in 2D TMy mode") + raise CmdInputError(f"'{self.params_str()}' polarisation must be y in 2D TMy mode") elif '2D TMz' in grid.mode and (polarisation == 'x' or polarisation == 'y'): - raise CmdInputError(f"'{self.__str__()}' polarisation must be z in 2D TMz mode") + raise CmdInputError(f"'{self.params_str()}' polarisation must be z in 2D TMz mode") - xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.__str__()) + xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.params_str()) # Check if there is a waveformID in the waveforms list if not any(x.ID == waveform_id for x in grid.waveforms): - raise CmdInputError(f"'{self.__str__()}' there is no waveform with the identifier {waveform_id}") + raise CmdInputError(f"'{self.params_str()}' there is no waveform with the identifier {waveform_id}") m = MagneticDipoleUser() m.polarisation = polarisation @@ -358,11 +358,11 @@ class MagneticDipole(UserObjectMulti): start = self.kwargs['start'] stop = self.kwargs['stop'] if start < 0: - raise CmdInputError(f"'{self.__str__()}' delay of the initiation of the source should not be less than zero") + raise CmdInputError(f"'{self.params_str()}' delay of the initiation of the source should not be less than zero") if stop < 0: - raise CmdInputError(f"'{self.__str__()}' time to remove the source should not be less than zero") + raise CmdInputError(f"'{self.params_str()}' time to remove the source should not be less than zero") if stop - start <= 0: - raise CmdInputError(f"'{self.__str__()}' duration of the source should not be zero or less") + raise CmdInputError(f"'{self.params_str()}' duration of the source should not be zero or less") m.start = start if stop > grid.timewindow: m.stop = grid.timewindow @@ -427,7 +427,7 @@ class TransmissionLine(UserObjectMulti): elif '2D TMz' in grid.mode and (polarisation == 'x' or polarisation == 'y'): raise CmdInputError(f"'{self.params_str()}' polarisation must be z in 2D TMz mode") - xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.__str__()) + xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.params_str()) if resistance <= 0 or resistance >= z0: raise CmdInputError(f"'{self.params_str()}' requires a resistance greater than zero and less than the impedance of free space, i.e. 376.73 Ohms") @@ -498,9 +498,9 @@ class Rx(UserObjectMulti): try: p1 = self.kwargs['p1'] except KeyError: - raise CmdInputError(f"'{self.__str__()}' has an incorrect number of parameters") + raise CmdInputError(f"'{self.params_str()}' has an incorrect number of parameters") - p = uip.check_src_rx_point(p1, self.__str__()) + p = uip.check_src_rx_point(p1, self.params_str()) r = self.constructor() r.xcoord, r.ycoord, r.zcoord = p @@ -519,7 +519,7 @@ class Rx(UserObjectMulti): if field in allowableoutputs: r.outputs[field] = np.zeros(grid.iterations, dtype=config.dtypes['float_or_double']) else: - raise CmdInputError(f"'{self.__str__()}' contains an output type that is not allowable. Allowable outputs in current context are {allowableoutputs}") + raise CmdInputError(f"'{self.params_str()}' contains an output type that is not allowable. Allowable outputs in current context are {allowableoutputs}") # If no ID or outputs are specified, use default except KeyError: r.ID = r.__class__.__name__ + '(' + str(r.xcoord) + ',' + str(r.ycoord) + ',' + str(r.zcoord) + ')' @@ -555,31 +555,31 @@ class RxArray(UserObjectMulti): p2 = self.kwargs['p2'] dl = self.kwargs['dl'] except KeyError: - raise CmdInputError(f"'{self.__str__()}' requires exactly 9 parameters") + raise CmdInputError(f"'{self.params_str()}' requires exactly 9 parameters") - xs, ys, zs = uip.check_src_rx_point(p1, self.__str__(), 'lower') - xf, yf, zf = uip.check_src_rx_point(p2, self.__str__(), 'upper') + xs, ys, zs = uip.check_src_rx_point(p1, self.params_str(), 'lower') + xf, yf, zf = uip.check_src_rx_point(p2, self.params_str(), 'upper') dx, dy, dz = uip.discretise_point(dl) if xs > xf or ys > yf or zs > zf: - raise CmdInputError(f"'{self.__str__()}' the lower coordinates should be less than the upper coordinates") + raise CmdInputError(f"'{self.params_str()}' the lower coordinates should be less than the upper coordinates") if dx < 0 or dy < 0 or dz < 0: - raise CmdInputError(f"'{self.__str__()}' the step size should not be less than zero") + raise CmdInputError(f"'{self.params_str()}' the step size should not be less than zero") if dx < 1: if dx == 0: dx = 1 else: - raise CmdInputError(f"'{self.__str__()}' the step size should not be less than the spatial discretisation") + raise CmdInputError(f"'{self.params_str()}' the step size should not be less than the spatial discretisation") if dy < 1: if dy == 0: dy = 1 else: - raise CmdInputError(f"'{self.__str__()}' the step size should not be less than the spatial discretisation") + raise CmdInputError(f"'{self.params_str()}' the step size should not be less than the spatial discretisation") if dz < 1: if dz == 0: dz = 1 else: - raise CmdInputError(f"'{self.__str__()}' the step size should not be less than the spatial discretisation") + raise CmdInputError(f"'{self.params_str()}' the step size should not be less than the spatial discretisation") log.info(f'Receiver array {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m with steps {dx * grid.dx:g}m, {dy * grid.dy:g}m, {dz * grid.dz:g}m') @@ -624,16 +624,16 @@ class Snapshot(UserObjectMulti): def create(self, grid, uip): if isinstance(grid, SubGridBase): - raise CmdInputError(f"'{self.__str__()}' do not add Snapshots to Subgrids.") + raise CmdInputError(f"'{self.params_str()}' do not add Snapshots to Subgrids.") try: p1 = self.kwargs['p1'] p2 = self.kwargs['p2'] dl = self.kwargs['dl'] filename = self.kwargs['filename'] except KeyError: - raise CmdInputError(f"'{self.__str__()}' requires exactly 11 parameters") + raise CmdInputError(f"'{self.params_str()}' requires exactly 11 parameters") - p1, p2 = uip.check_box_points(p1, p2, self.__str__()) + p1, p2 = uip.check_box_points(p1, p2, self.params_str()) xs, ys, zs = p1 xf, yf, zf = p2 dx, dy, dz = uip.discretise_point(dl) @@ -646,18 +646,18 @@ class Snapshot(UserObjectMulti): try: time = self.kwargs['time'] except KeyError: - raise CmdInputError(f"'{self.__str__()}' requires exactly 5 parameters") + raise CmdInputError(f"'{self.params_str()}' requires exactly 5 parameters") if time > 0: iterations = round_value((time / grid.dt)) + 1 else: - raise CmdInputError(f"'{self.__str__()}' time value must be greater than zero") + raise CmdInputError(f"'{self.params_str()}' time value must be greater than zero") if dx < 0 or dy < 0 or dz < 0: - raise CmdInputError(f"'{self.__str__()}' the step size should not be less than zero") + raise CmdInputError(f"'{self.params_str()}' the step size should not be less than zero") if dx < 1 or dy < 1 or dz < 1: - raise CmdInputError(f"'{self.__str__()}' the step size should not be less than the spatial discretisation") + raise CmdInputError(f"'{self.params_str()}' the step size should not be less than the spatial discretisation") if iterations <= 0 or iterations > grid.iterations: - raise CmdInputError(f"'{self.__str__()}' time value is not valid") + raise CmdInputError(f"'{self.params_str()}' time value is not valid") # Replace with old style snapshots if there are subgrids #if grid.subgrids: @@ -1005,7 +1005,7 @@ class GeometryView(UserObjectMulti): GeometryViewUser = self.geometry_view_constructor(grid, output_type) - p1, p2 = uip.check_box_points(p1, p2, self.__str__()) + p1, p2 = uip.check_box_points(p1, p2, self.params_str()) xs, ys, zs = p1 xf, yf, zf = p2 @@ -1061,7 +1061,7 @@ class GeometryObjectsWrite(UserObjectMulti): except KeyError: raise CmdInputError(f"'{self.params_str()}' requires exactly seven parameters") - p1, p2 = uip.check_box_points(p1, p2, self.__str__()) + p1, p2 = uip.check_box_points(p1, p2, self.params_str()) x0, y0, z0 = p1 x1, y1, z1 = p2 diff --git a/gprMax/cmds_single_use.py b/gprMax/cmds_single_use.py index f92b8762..e6ee764b 100644 --- a/gprMax/cmds_single_use.py +++ b/gprMax/cmds_single_use.py @@ -34,7 +34,7 @@ import gprMax.config as config from .exceptions import CmdInputError from .waveforms import Waveform from .utilities import round_value - +from .utilities import set_openmp_threads log = logging.getLogger(__name__) @@ -78,11 +78,9 @@ class Domain(UserObjectSingle): def __str__(self): try: - s = '#domain: {} {} {}'.format(self.kwargs['p1'][0], - self.kwargs['p1'][1], - self.kwargs['p1'][2]) + s = f"#domain: {self.kwargs['p1'][0]} {self.kwargs['p1'][1]} {self.kwargs['p1'][2]}" except KeyError: - log.warning('error message') + log.exception('error message') return s @@ -95,7 +93,7 @@ class Domain(UserObjectSingle): if G.nx == 0 or G.ny == 0 or G.nz == 0: raise CmdInputError(f"'{self.params_str()}' requires at least one cell in every dimension") - log.info(f'Domain size: {self.kwargs['p1'][0]:g} x {self.kwargs['p1'][1]:g} x {self.kwargs['p1'][2]:g}m ({G.nx:d} x {G.ny:d} x {G.nz:d} = {(G.nx * G.ny * G.nz):g} cells)') + log.info(f"Domain size: {self.kwargs['p1'][0]:g} x {self.kwargs['p1'][1]:g} x {self.kwargs['p1'][2]:g}m ({G.nx:d} x {G.ny:d} x {G.nz:d} = {(G.nx * G.ny * G.nz):g} cells)") # Calculate time step at CFL limit; switch off appropriate PMLs for 2D if G.nx == 1: @@ -121,33 +119,18 @@ class Domain(UserObjectSingle): log.info(f'Time step (at CFL limit): {G.dt:g} secs') # Number of threads (OpenMP) to use - if sys.platform == 'darwin': - os.environ['OMP_WAIT_POLICY'] = 'ACTIVE' # Should waiting threads consume CPU power (can drastically effect performance) - os.environ['OMP_DYNAMIC'] = 'FALSE' # Number of threads may be adjusted by the run time environment to best utilize system resources - os.environ['OMP_PLACES'] = 'cores' # Each place corresponds to a single core (having one or more hardware threads) - os.environ['OMP_PROC_BIND'] = 'TRUE' # Bind threads to physical cores - # os.environ['OMP_DISPLAY_ENV'] = 'TRUE' # Prints OMP version and environment variables (useful for debug) - - # Catch bug with Windows Subsystem for Linux (https://github.com/Microsoft/BashOnWindows/issues/785) - if 'Microsoft' in config.hostinfo['osversion']: - os.environ['KMP_AFFINITY'] = 'disabled' - del os.environ['OMP_PLACES'] - del os.environ['OMP_PROC_BIND'] - - if os.environ.get('OMP_NUM_THREADS'): - G.nthreads = int(os.environ.get('OMP_NUM_THREADS')) - else: - # Set number of threads to number of physical CPU cores - G.nthreads = config.hostinfo['physicalcores'] - os.environ['OMP_NUM_THREADS'] = str(G.nthreads) - + G.nthreads = set_openmp_threads() log.info(f'Number of CPU (OpenMP) threads: {G.nthreads}') + if G.nthreads > config.hostinfo['physicalcores']: - log.warning(Fore.RED + f'You have specified more threads ({G.nthreads}) than available physical CPU cores ({config.hostinfo['physicalcores']}). This may lead to degraded performance.' + Style.RESET_ALL) + log.warning(Fore.RED + f"You have specified more threads ({G.nthreads}) \ + than available physical CPU cores ({config.hostinfo['physicalcores']}). \ + This may lead to degraded performance." + Style.RESET_ALL) class Discretisation(UserObjectSingle): - """Allows you to specify the discretization of space in the x , y and z directions respectively + """Allows you to specify the discretization of space in the x, y, and z + directions respectively. :param p1: Specify discretisation in x, y, z direction :type p1: list of floats, non-optional @@ -159,11 +142,9 @@ class Discretisation(UserObjectSingle): def __str__(self): try: - s = '#dx_dy_dz: {} {} {}'.format(self.kwargs['p1'][0], - self.kwargs['p1'][1], - self.kwargs['p1'][2]) + s = f"#dx_dy_dz: {self.kwargs['p1'][0]} {self.kwargs['p1'][1]} {self.kwargs['p1'][2]}" except KeyError: - log.info('error message') + log.exception('error message') return s @@ -172,17 +153,16 @@ class Discretisation(UserObjectSingle): G.dl = np.array(self.kwargs['p1']) G.dx, G.dy, G.dz = self.kwargs['p1'] except KeyError: - raise CmdInputError('Discretisation requires a point') + raise CmdInputError(f"'{self.params_str()}' discretisation requires a point") if G.dl[0] <= 0: - raise CmdInputError('Discretisation requires the x-direction spatial step to be greater than zero') + raise CmdInputError(f"'{self.params_str()}' discretisation requires the x-direction spatial step to be greater than zero") if G.dl[1] <= 0: - raise CmdInputError(' Discretisation requires the y-direction spatial step to be greater than zero') + raise CmdInputError(f"'{self.params_str()}' discretisation requires the y-direction spatial step to be greater than zero") if G.dl[2] <= 0: - raise CmdInputError('Discretisation requires the z-direction spatial step to be greater than zero') + raise CmdInputError(f"'{self.params_str()}' discretisation requires the z-direction spatial step to be greater than zero") - if config.is_messages(): - log.info('Spatial discretisation: {:g} x {:g} x {:g}m'.format(*G.dl)) + log.info(f'Spatial discretisation: {G.dl[0]:g} x {G.dl[1]:g} x {G.dl[2]:g}m') class TimeWindow(UserObjectSingle): @@ -200,12 +180,12 @@ class TimeWindow(UserObjectSingle): def __str__(self): try: - s = '#time_window: {}'.format(self.kwargs['time']) + s = f"#time_window: {self.kwargs['time']}" except KeyError: try: - s = '#time_window: {}'.format(self.kwargs['iterations']) + s = f"#time_window: {self.kwargs['iterations']}" except KeyError: - log.info('time window error') + log.exception('time window error') return s @@ -217,7 +197,6 @@ class TimeWindow(UserObjectSingle): iterations = int(self.kwargs['iterations']) G.timewindow = (iterations - 1) * G.dt G.iterations = iterations - except KeyError: pass @@ -227,16 +206,14 @@ class TimeWindow(UserObjectSingle): G.timewindow = tmp G.iterations = int(np.ceil(tmp / G.dt)) + 1 else: - raise CmdInputError(self.__str__() + ' must have a value greater than zero') - + raise CmdInputError(f"'{self.params_str()}' must have a value greater than zero") except KeyError: pass if not G.timewindow: - raise CmdInputError('TimeWindow: Specify a time or number of iterations') + raise CmdInputError(f"'{self.params_str()}' specify a time or number of iterations") - if config.is_messages(): - log.info('Time window: {:g} secs ({} iterations)'.format(G.timewindow, G.iterations)) + log.info(f'Time window: {G.timewindow:g} secs ({G.iterations} iterations)') class Messages(UserObjectSingle): diff --git a/gprMax/contexts.py b/gprMax/contexts.py index 7ca31ba0..6a902689 100644 --- a/gprMax/contexts.py +++ b/gprMax/contexts.py @@ -70,8 +70,8 @@ class NoMPIContext(Context): is parallelised using either OpenMP (CPU) or CUDA (GPU). """ - def __init__(self): - super().__init__() + def __init__(self, sim_config): + super().__init__(sim_config) def _run(self): """Specialise how the models are farmed out.""" @@ -96,13 +96,13 @@ class NoMPIContext(Context): model.solve(solver) def make_time_report(self): - """Function to specialise the time reporting for the standard Simulation + """Function to specialise the time reporting for the standard simulation context. """ sim_time = datetime.timedelta(seconds=self.tsimend - self.tsimstart) - s = f'\n=== Simulation on {self.simconfig.hostinfo['hostname']} completed in [HH:MM:SS]: {sim_time}' - return f'{s} {'=' * (get_terminal_width() - 1 - len(s))}\n' + s = f"\n=== Simulation on {self.simconfig.hostinfo['hostname']} completed in [HH:MM:SS]: {sim_time}" + return f"{s} {'=' * (get_terminal_width() - 1 - len(s))}\n" class MPIContext(Context): @@ -111,9 +111,8 @@ class MPIContext(Context): or CUDA (GPU). """ - def __init__(self): - super().__init__() - + def __init__(self, sim_config): + super().__init__(sim_config) from mpi4py import MPI def _run(self): @@ -125,8 +124,9 @@ class MPIContext(Context): class MPINoSpawnContext(Context): - def __init__(self): - super().__init__() + def __init__(self, sim_config): + super().__init__(sim_config) + from mpi4py import MPI def _run(self): pass diff --git a/gprMax/exceptions.py b/gprMax/exceptions.py index 1beab97b..57a1c15b 100644 --- a/gprMax/exceptions.py +++ b/gprMax/exceptions.py @@ -16,13 +16,14 @@ # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . +import logging import sys from colorama import init from colorama import Fore - init() +log = logging.getLogger(__name__) sys.tracebacklimit = None @@ -33,7 +34,7 @@ class GeneralError(ValueError): self.message = message super(GeneralError, self).__init__(message, *args) - print(Fore.RED) + log.exception(Fore.RED) class CmdInputError(ValueError): @@ -45,4 +46,4 @@ class CmdInputError(ValueError): self.message = message super(CmdInputError, self).__init__(message, *args) - print(Fore.RED) + log.exception(Fore.RED) diff --git a/gprMax/fields_outputs.py b/gprMax/fields_outputs.py index decaab96..87d5cb4d 100644 --- a/gprMax/fields_outputs.py +++ b/gprMax/fields_outputs.py @@ -31,7 +31,7 @@ def store_outputs(G): iteration (int): Current iteration number. Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic field values. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ iteration = G.iteration @@ -105,7 +105,7 @@ def write_hdf5_main_grid_outputfile(outputfile, G): Args: outputfile (str): Name of the output file. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ write_data(outputfile, G) @@ -116,7 +116,7 @@ def write_hdf5_sub_grid_outputfile(outputfile, G): Args: outputfile (str): Name of the output file. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ stem = outputfile.stem @@ -145,7 +145,7 @@ def write_data(outputfile, G): Args: outputfile (str): Name of the output file. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ f = h5py.File(outputfile, 'w') diff --git a/gprMax/fractals.py b/gprMax/fractals.py index f2ff73c2..7352e056 100644 --- a/gprMax/fractals.py +++ b/gprMax/fractals.py @@ -27,7 +27,7 @@ from .utilities import round_value np.seterr(divide='raise') -class FractalSurface(object): +class FractalSurface: """Fractal surfaces.""" surfaceIDs = ['xminus', 'xplus', 'yminus', 'yplus', 'zminus', 'zplus'] @@ -109,7 +109,7 @@ class FractalSurface(object): - self.fractalrange[0]) / fractalrange) * fractalmin) -class FractalVolume(object): +class FractalVolume: """Fractal volumes.""" def __init__(self, xs, xf, ys, yf, zs, zf, dimension): @@ -209,7 +209,7 @@ class FractalVolume(object): self.mask[maskxs:maskxf, maskys:maskyf, maskzs:maskzf] = 1 -class Grass(object): +class Grass: """Geometry information for blades of grass.""" def __init__(self, numblades): diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py index 78ddf4e5..56321c9b 100644 --- a/gprMax/geometry_outputs.py +++ b/gprMax/geometry_outputs.py @@ -31,7 +31,7 @@ from .cython.geometry_outputs import define_fine_geometry from .utilities import round_value -class GeometryView(object): +class GeometryView: """Views of the geometry of the model.""" if sys.byteorder == 'little': @@ -313,7 +313,7 @@ class GeometryView(object): f.write('\n'.encode('utf-8')) -class GeometryObjects(object): +class GeometryObjects: """Geometry objects to be written to file.""" def __init__(self, xs=None, ys=None, zs=None, xf=None, yf=None, zf=None, basefilename=None): @@ -390,12 +390,13 @@ class GeometryObjects(object): fmaterials.write(dispersionstr + '\n') class GeometryViewFineMultiGrid: - """Geometry view for all grids in the simulation.""" + """Geometry view for all grids in the simulation. - """"Slicing is not supported by this class :( - only the full extent of the grids - are output. The subgrids are output without the non-working regions If you - require domainslicing GeometryView seperately for each grid you require and - view them at once in Paraview.""" + Slicing is not supported by this class :( - only the full extent of the grids + are output. The subgrids are output without the non-working regions If you + require domainslicing GeometryView seperately for each grid you require and + view them at once in Paraview. + """ if sys.byteorder == 'little': byteorder = 'LittleEndian' diff --git a/gprMax/grid.py b/gprMax/grid.py index d6b9f784..b24bf3a2 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -24,8 +24,6 @@ from colorama import Fore from colorama import Style init() import numpy as np -np.seterr(invalid='raise') -from scipy.constants import c import gprMax.config as config from .exceptions import GeneralError @@ -35,6 +33,8 @@ from .utilities import fft_power from .utilities import human_size from .utilities import round_value +np.seterr(invalid='raise') + class FDTDGrid: """Holds attributes associated with entire grid. A convenient way for @@ -211,7 +211,8 @@ class FDTDGrid: and/or run model on the host. """ if self.memoryusage > config.hostinfo['ram']: - raise GeneralError(f'Memory (RAM) required ~{human_size(self.memoryusage)} exceeds {human_size(config.hostinfo['ram'], a_kilobyte_is_1024_bytes=True)} detected!\n') + raise GeneralError(f"Memory (RAM) required ~{human_size(self.memoryusage)} \ + exceeds {human_size(config.hostinfo['ram'], a_kilobyte_is_1024_bytes=True)} detected!\n") def tmx(self): """Add PEC boundaries to invariant direction in 2D TMx mode. @@ -269,9 +270,14 @@ class CUDAGrid(FDTDGrid): """Additional grid methods for solving on GPU using CUDA.""" def set_blocks_per_grid(self): - """Set the blocks per grid size used for updating the electric and magnetic field arrays on a GPU.""" + """Set the blocks per grid size used for updating the electric and + magnetic field arrays on a GPU. + """ - config.cuda['gpus'].bpg = (int(np.ceil(((self.nx + 1) * (self.ny + 1) * (self.nz + 1)) / config.cuda['gpus'].tpb[0])), 1, 1) + config.cuda['gpus'].bpg = (int(np.ceil(((self.nx + 1) * + (self.ny + 1) * + (self.nz + 1)) / + config.cuda['gpus'].tpb[0])), 1, 1) def initialise_arrays(self): """Initialise geometry and field arrays on GPU.""" @@ -303,7 +309,9 @@ class CUDAGrid(FDTDGrid): if config.cuda['gpus'] is not None: if self.memoryusage - snapsmemsize > config.cuda['gpus'].totalmem: - raise GeneralError(f'Memory (RAM) required ~{human_size(self.memoryusage)} exceeds {human_size(config.cuda['gpus'].totalmem, a_kilobyte_is_1024_bytes=True} detected on specified {config.cuda['gpus'].deviceID} - {config.cuda['gpus'].name} GPU!\n') + raise GeneralError(f"Memory (RAM) required ~{human_size(self.memoryusage)} \ + exceeds {human_size(config.cuda['gpus'].totalmem, a_kilobyte_is_1024_bytes=True)} \ + detected on specified {config.cuda['gpus'].deviceID} - {config.cuda['gpus'].name} GPU!\n") # If the required memory for the model without the snapshots will # fit on the GPU then transfer and store snaphots on host @@ -439,7 +447,7 @@ def Ix(x, y, z, Hx, Hy, Hz, G): Args: x, y, z (float): Coordinates of position in grid. Hx, Hy, Hz (memory view): numpy array of magnetic field values. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ if y == 0 or z == 0: @@ -456,7 +464,7 @@ def Iy(x, y, z, Hx, Hy, Hz, G): Args: x, y, z (float): Coordinates of position in grid. Hx, Hy, Hz (memory view): numpy array of magnetic field values. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ if x == 0 or z == 0: @@ -473,7 +481,7 @@ def Iz(x, y, z, Hx, Hy, Hz, G): Args: x, y, z (float): Coordinates of position in grid. Hx, Hy, Hz (memory view): numpy array of magnetic field values. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ if x == 0 or y == 0: diff --git a/gprMax/input_cmds_file.py b/gprMax/hash_cmds_file.py similarity index 93% rename from gprMax/input_cmds_file.py rename to gprMax/hash_cmds_file.py index 89bcef9a..5a20f273 100644 --- a/gprMax/input_cmds_file.py +++ b/gprMax/hash_cmds_file.py @@ -16,15 +16,18 @@ # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . +from io import StringIO +import logging import os import sys -from io import StringIO import gprMax.config as config from .exceptions import CmdInputError -from .input_cmds_geometry import process_geometrycmds -from .input_cmds_multiuse import process_multicmds -from .input_cmds_singleuse import process_singlecmds +from .hash_cmds_geometry import process_geometrycmds +from .hash_cmds_multiuse import process_multicmds +from .hash_cmds_singleuse import process_singlecmds + +log = logging.getLogger(__name__) def process_python_include_code(inputfile, usernamespace): @@ -95,7 +98,7 @@ def process_python_include_code(inputfile, usernamespace): # Print any generated output that is not commands if pythonout: - print('Python messages (from stdout/stderr): {}\n'.format(pythonout)) + log.info(f'Python messages (from stdout/stderr): {pythonout}\n') # Add any other commands to list elif(inputlines[x].startswith('#')): @@ -164,13 +167,15 @@ def write_processed_file(processedlines, appendmodelnumber): appendmodelnumber (str): Text to append to filename. """ - processedfile = os.path.join(config.general['outputfilepath'], os.path.splitext(config.general['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: f.write('{}'.format(item)) - print('Written input commands, after processing any Python code and include commands, to file: {}\n'.format(processedfile)) + log.info(f'Written input commands, after processing any Python code and include commands, to file: {processedfile}\n') def check_cmd_names(processedlines, checkessential=True): @@ -192,7 +197,7 @@ def check_cmd_names(processedlines, checkessential=True): essentialcmds = ['#domain', '#dx_dy_dz', '#time_window'] # Commands that there should only be one instance of in a model - singlecmds = dict.fromkeys(['#domain', '#dx_dy_dz', '#time_window', '#title', '#messages', '#num_threads', '#time_step_stability_factor', '#pml_formulation', '#pml_cells', '#excitation_file', '#src_steps', '#rx_steps', '#taguchi', '#end_taguchi', '#output_dir'], None) + singlecmds = dict.fromkeys(['#domain', '#dx_dy_dz', '#time_window', '#title', '#messages', '#num_threads', '#time_step_stability_factor', '#pml_formulation', '#pml_cells', '#excitation_file', '#src_steps', '#rx_steps', '#output_dir'], None) # Commands that there can be multiple instances of in a model - these will be lists within the dictionary multiplecmds = {key: [] for key in ['#geometry_view', '#geometry_objects_write', '#material', '#soil_peplinski', '#add_dispersion_debye', '#add_dispersion_lorentz', '#add_dispersion_drude', '#waveform', '#voltage_source', '#hertzian_dipole', '#magnetic_dipole', '#transmission_line', '#rx', '#rx_array', '#snapshot', '#pml_cfs', '#include_file']} @@ -282,7 +287,7 @@ def parse_hash_commands(model_config, G, scene): Args: model_config (ModelConfig): Model level configuration object. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. scene (Scene): Scene object. Returns: @@ -302,9 +307,8 @@ def parse_hash_commands(model_config, G, scene): uservars = '' for key, value in sorted(usernamespace.items()): if key != '__builtins__': - uservars += '{}: {}, '.format(key, value) - usv_s = """Constants/variables used/available for Python scripting: {{{}}}\n""" - print(usv_s.format(uservars[:-2])) + uservars += f'{key}: {value}, ' + log.info(f'Constants/variables used/available for Python scripting: {{{uservars[:-2]}}}\n') # Write a file containing the input commands after Python or include # file commands have been processed diff --git a/gprMax/input_cmds_geometry.py b/gprMax/hash_cmds_geometry.py similarity index 96% rename from gprMax/input_cmds_geometry.py rename to gprMax/hash_cmds_geometry.py index 74e8ca36..800ceede 100644 --- a/gprMax/input_cmds_geometry.py +++ b/gprMax/hash_cmds_geometry.py @@ -16,7 +16,6 @@ # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -import sys import numpy as np from .cmds_geometry.geometry_objects_read import GeometryObjectsRead @@ -233,13 +232,11 @@ def process_geometrycmds(geometry): n_materials = round_value(tmp[11]) mixing_model_id = tmp[12] ID = tmp[13] - # without seed + if len(tmp) == 14: fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, id=ID, n_materials=n_materials) - # with seed elif len(tmp) == 15: fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, id=ID, n_materials=n_materials, seed=tmp[14]) - # user specified averaging elif len(tmp) == 16: fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, id=ID, n_materials=n_materials, seed=tmp[14], averaging=tmp[15].lower()) else: @@ -262,7 +259,6 @@ def process_geometrycmds(geometry): limits = [float(tmp[10]), float(tmp[11])] fractal_box_id = tmp[12] - # No seed if len(tmp) == 13: asr = AddSurfaceRoughness(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, limits=limits, fractal_box_id=fractal_box_id) elif len(tmp) == 14: @@ -295,13 +291,13 @@ def process_geometrycmds(geometry): n_blades = int(tmp[10]) fractal_box_id = tmp[11] - # no seed if len(tmp) == 12: grass = AddGrass(p1=p1, p2=p2, frac_dim=frac_dim, limits=limits, n_blades=n_blades, fractal_box_id=fractal_box_id) elif len(tmp) == 13: grass = AddGrass(p1=p1, p2=p2, frac_dim=frac_dim, limits=limits, n_blades=n_blades, fractal_box_id=fractal_box_id, seed=int(tmp[12])) else: raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given') + scene_objects.append(grass) return scene_objects diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/hash_cmds_multiuse.py similarity index 90% rename from gprMax/input_cmds_multiuse.py rename to gprMax/hash_cmds_multiuse.py index 7098e672..527cc9e7 100644 --- a/gprMax/input_cmds_multiuse.py +++ b/gprMax/hash_cmds_multiuse.py @@ -15,23 +15,6 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -# Copyright (C) 2015-2018: The University of Edinburgh -# Authors: Craig Warren and Antonis Giannopoulos -# -# This file is part of gprMax. -# -# gprMax is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# gprMax is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with gprMax. If not, see . from .cmds_multiple import Waveform from .cmds_multiple import VoltageSource @@ -65,7 +48,6 @@ def process_multicmds(multicmds): scene_objects = [] - # Waveform definitions cmdname = '#waveform' if multicmds[cmdname] is not None: for cmdinstance in multicmds[cmdname]: @@ -76,7 +58,6 @@ def process_multicmds(multicmds): waveform = Waveform(wave_type=tmp[0], amp=float(tmp[1]), freq=float(tmp[2]), id=tmp[3]) scene_objects.append(waveform) - # Voltage source cmdname = '#voltage_source' if multicmds[cmdname] is not None: for cmdinstance in multicmds[cmdname]: @@ -90,7 +71,6 @@ def process_multicmds(multicmds): scene_objects.append(voltage_source) - # Hertzian dipole cmdname = '#hertzian_dipole' if multicmds[cmdname] is not None: for cmdinstance in multicmds[cmdname]: @@ -106,7 +86,6 @@ def process_multicmds(multicmds): scene_objects.append(hertzian_dipole) - # Magnetic dipole cmdname = '#magnetic_dipole' if multicmds[cmdname] is not None: for cmdinstance in multicmds[cmdname]: @@ -122,7 +101,6 @@ def process_multicmds(multicmds): scene_objects.append(magnetic_dipole) - # Transmission line cmdname = '#transmission_line' if multicmds[cmdname] is not None: for cmdinstance in multicmds[cmdname]: @@ -139,7 +117,6 @@ def process_multicmds(multicmds): scene_objects.append(tl) - # Receiver cmdname = '#rx' if multicmds[cmdname] is not None: for cmdinstance in multicmds[cmdname]: @@ -150,9 +127,9 @@ def process_multicmds(multicmds): rx = Rx(p1=(float(tmp[0]), float(tmp[1]), float(tmp[2]))) else: rx = Rx(p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])), id=tmp[3], outputs=tmp[4:]) + scene_objects.append(rx) - # Receiver array cmdname = '#rx_array' if multicmds[cmdname] is not None: for cmdinstance in multicmds[cmdname]: @@ -167,7 +144,6 @@ def process_multicmds(multicmds): rx_array = RxArray(p1=p1, p2=p2, dl=dl) scene_objects.append(rx_array) - # Snapshot cmdname = '#snapshot' if multicmds[cmdname] is not None: for cmdinstance in multicmds[cmdname]: @@ -190,7 +166,6 @@ def process_multicmds(multicmds): scene_objects.append(snapshot) - # Materials cmdname = '#material' if multicmds[cmdname] is not None: for cmdinstance in multicmds[cmdname]: @@ -279,7 +254,6 @@ def process_multicmds(multicmds): id=tmp[6]) scene_objects.append(soil) - # Geometry views (creates VTK-based geometry files) cmdname = '#geometry_view' if multicmds[cmdname] is not None: for cmdinstance in multicmds[cmdname]: @@ -292,10 +266,8 @@ def process_multicmds(multicmds): dl = float(tmp[6]), float(tmp[7]), float(tmp[8]) geometry_view = GeometryView(p1=p1, p2=p2, dl=dl, filename=tmp[9], output_type=tmp[10]) - scene_objects.append(geometry_view) - # Geometry object(s) output cmdname = '#geometry_objects_write' if multicmds[cmdname] is not None: for cmdinstance in multicmds[cmdname]: @@ -308,8 +280,6 @@ def process_multicmds(multicmds): gow = GeometryObjectsWrite(p1=p1, p2=p2, filename=tmp[6]) scene_objects.append(gow) - - # Complex frequency shifted (CFS) PML parameter cmdname = '#pml_cfs' if multicmds[cmdname] is not None: if len(multicmds[cmdname]) > 2: @@ -335,8 +305,3 @@ def process_multicmds(multicmds): scene_objects.append(pml_cfs) return scene_objects - - -def process_subgrid_hsg(cmdinstance): - - pass diff --git a/gprMax/input_cmds_singleuse.py b/gprMax/hash_cmds_singleuse.py similarity index 85% rename from gprMax/input_cmds_singleuse.py rename to gprMax/hash_cmds_singleuse.py index a1efa76a..b0284c7f 100644 --- a/gprMax/input_cmds_singleuse.py +++ b/gprMax/hash_cmds_singleuse.py @@ -15,23 +15,6 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -# Copyright (C) 2015-2018: The University of Edinburgh -# Authors: Craig Warren and Antonis Giannopoulos -# -# This file is part of gprMax. -# -# gprMax is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# gprMax is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with gprMax. If not, see . from .cmds_single_use import Messages from .cmds_single_use import Title @@ -55,7 +38,6 @@ def process_singlecmds(singlecmds): Args: singlecmds (dict): Commands that can only occur once in the model. - G (Grid): Holds essential parameters describing the model. Returns: scene_objects (list): Holds objects in scene. @@ -78,7 +60,6 @@ def process_singlecmds(singlecmds): title = Title(name=str(singlecmds[cmd])) scene_objects.append(title) - # Set the output directory cmd = '#output_dir' if singlecmds[cmd] is not None: output_dir = OutputDir(dir=singlecmds[cmd]) @@ -94,7 +75,6 @@ def process_singlecmds(singlecmds): num_thread = NumThreads(n=tmp[0]) scene_objects.append(num_thread) - # Spatial discretisation cmd = '#dx_dy_dz' if singlecmds[cmd] is not None: tmp = [float(x) for x in singlecmds[cmd].split()] @@ -129,8 +109,6 @@ def process_singlecmds(singlecmds): tmp = tmp[0].lower() # If number of iterations given - # The +/- 1 used in calculating the number of iterations is to account for - # the fact that the solver (iterations) loop runs from 0 to < G.iterations try: tmp = int(tmp) tw = TimeWindow(iterations=tmp) @@ -163,6 +141,7 @@ def process_singlecmds(singlecmds): tmp = singlecmds[cmd].split() if len(tmp) != 3: raise CmdInputError(cmd + ' requires exactly three parameters') + p1 = (float(tmp[0]), float(tmp[1]), float(tmp[2])) src_steps = SrcSteps(p1=p1) scene_objects.append(src_steps) diff --git a/gprMax/materials.py b/gprMax/materials.py index d75f22b1..70e14f52 100644 --- a/gprMax/materials.py +++ b/gprMax/materials.py @@ -21,7 +21,7 @@ import numpy as np import gprMax.config as config -class Material(object): +class Material: """Super-class to describe generic, non-dispersive materials, their properties and update coefficients. """ @@ -49,7 +49,7 @@ class Material(object): """Calculates the magnetic update coefficients of the material. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ HA = (config.m0 * self.mr / G.dt) + 0.5 * self.sm @@ -64,7 +64,7 @@ class Material(object): """Calculates the electric update coefficients of the material. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ EA = (config.e0 * self.er / G.dt) + 0.5 * self.se @@ -126,7 +126,7 @@ class DispersiveMaterial(Material): """Calculates the electric update coefficients of the material. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ # The implementation of the dispersive material modelling comes from the @@ -208,7 +208,7 @@ def process_materials(G): store in arrays, and build text list of materials/properties Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. Returns: materialsdata (list): List of material IDs, names, and properties to @@ -274,7 +274,7 @@ def process_materials(G): return materialsdata -class PeplinskiSoil(object): +class PeplinskiSoil: """Soil objects that are characterised according to a mixing model by Peplinski (http://dx.doi.org/10.1109/36.387598). """ @@ -307,7 +307,7 @@ class PeplinskiSoil(object): Args: nbins (int): Number of bins to use to create the different materials. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ # Debye model properties of water @@ -374,7 +374,7 @@ def create_built_in_materials(G): """Create pre-defined (built-in) materials. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ G.n_built_in_materials = len(G.materials) diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py index 7545af3a..f713801e 100644 --- a/gprMax/model_build_run.py +++ b/gprMax/model_build_run.py @@ -45,13 +45,13 @@ from .fields_outputs import kernel_template_store_outputs from .fields_outputs import write_hdf5_outputfile from .grid import FDTDGrid from .grid import dispersion_analysis -from .input_cmds_geometry import process_geometrycmds -from .input_cmds_file import process_python_include_code -from .input_cmds_file import write_processed_file -from .input_cmds_file import check_cmd_names -from .input_cmds_file import parse_hash_commands -from .input_cmds_singleuse import process_singlecmds -from .input_cmds_multiuse import process_multicmds +from .hash_cmds_geometry import process_geometrycmds +from .hash_cmds_file import process_python_include_code +from .hash_cmds_file import write_processed_file +from .hash_cmds_file import check_cmd_names +from .hash_cmds_file import parse_hash_commands +from .hash_cmds_singleuse import process_singlecmds +from .hash_cmds_multiuse import process_multicmds from .materials import Material from .materials import process_materials from .pml import CFS @@ -72,11 +72,10 @@ from .utilities import human_size from .utilities import open_path_file from .utilities import round32 from .utilities import timer -from .utilities import Printer - log = logging.getLogger(__name__) + class ModelBuildRun: """Builds and runs (solves) a model.""" @@ -84,7 +83,6 @@ class ModelBuildRun: self.G = G self.sim_config = sim_config self.model_config = model_config - self.printer = Printer(config) # Monitor memory usage self.p = None @@ -122,7 +120,7 @@ class ModelBuildRun: # Write files for any geometry views and geometry object outputs if not (G.geometryviews or G.geometryobjectswrite) and self.sim_config.geometry_only and config.is_messages(): - print(Fore.RED + '\nWARNING: No geometry views or geometry objects to output found.' + Style.RESET_ALL) + log.warning(Fore.RED + f'\nNo geometry views or geometry objects to output found.' + Style.RESET_ALL) if config.is_messages(): print() for i, geometryview in enumerate(G.geometryviews): geometryview.set_filename(self.model_config.appendmodelnumber) @@ -144,8 +142,7 @@ class ModelBuildRun: sim_config = self.sim_config G = self.G - printer = Printer(config) - printer.print(model_config.next_model) + log.info(model_config.next_model) scene = self.build_scene() @@ -172,9 +169,11 @@ class ModelBuildRun: config.materials['dispersiveCdtype'] = config.dtypes['C_float_or_double'] # Update estimated memory (RAM) usage - G.memoryusage += int(3 * config.materials['maxpoles'] * (G.nx + 1) * (G.ny + 1) * (G.nz + 1) * np.dtype(config.materials['dispersivedtype']).itemsize) + G.memoryusage += int(3 * config.materials['maxpoles'] * + (G.nx + 1) * (G.ny + 1) * (G.nz + 1) * + np.dtype(config.materials['dispersivedtype']).itemsize) G.memory_check() - printer.print('\nMemory (RAM) required - updated (dispersive): ~{}\n'.format(human_size(G.memoryusage))) + log.info(f'\nMemory (RAM) required - updated (dispersive): ~{human_size(G.memoryusage)}\n') for gb in gridbuilders: gb.grid.initialise_dispersive_arrays(config.materials['dispersivedtype']) @@ -188,7 +187,7 @@ class ModelBuildRun: G.memoryusage += int(snapsmemsize) G.memory_check(snapsmemsize=int(snapsmemsize)) - printer.print('\nMemory (RAM) required - updated (snapshots): ~{}\n'.format(human_size(G.memoryusage))) + log.info(f'\nMemory (RAM) required - updated (snapshots): ~{human_size(G.memoryusage)}\n') # Build materials for gb in gridbuilders: @@ -197,20 +196,33 @@ class ModelBuildRun: # Check to see if numerical dispersion might be a problem results = dispersion_analysis(G) if results['error']: - printer.print(Fore.RED + "\nWARNING: Numerical dispersion analysis not carried out as {}".format(results['error']) + Style.RESET_ALL) + log.warning(Fore.RED + f"\nNumerical dispersion analysis not carried \ + out as {results['error']}" + Style.RESET_ALL) elif results['N'] < config.numdispersion['mingridsampling']: - raise GeneralError("Non-physical wave propagation: Material '{}' has wavelength sampled by {} cells, less than required minimum for physical wave propagation. Maximum significant frequency estimated as {:g}Hz".format(results['material'].ID, results['N'], results['maxfreq'])) + raise GeneralError(f"Non-physical wave propagation: Material \ + '{results['material'].ID}' has wavelength sampled \ + by {results['N']} cells, less than required \ + minimum for physical wave propagation. Maximum \ + significant frequency estimated as {results['maxfreq']:g}Hz") elif results['deltavp'] and np.abs(results['deltavp']) > config.numdispersion['maxnumericaldisp']: - printer.print(Fore.RED + "\nWARNING: Potentially significant numerical dispersion. Estimated largest physical phase-velocity error is {:.2f}% in material '{}' whose wavelength sampled by {} cells. Maximum significant frequency estimated as {:g}Hz".format(results['deltavp'], results['material'].ID, results['N'], results['maxfreq']) + Style.RESET_ALL) + log.warning(Fore.RED + f"\nPotentially significant numerical dispersion. \ + Estimated largest physical phase-velocity error is \ + {results['deltavp']:.2f}% in material '{results['material'].ID}' \ + whose wavelength sampled by {results['N']} cells. Maximum \ + significant frequency estimated as {results['maxfreq']:g}Hz" + Style.RESET_ALL) elif results['deltavp']: - printer.print("\nNumerical dispersion analysis: estimated largest physical phase-velocity error is {:.2f}% in material '{}' whose wavelength sampled by {} cells. Maximum significant frequency estimated as {:g}Hz".format(results['deltavp'], results['material'].ID, results['N'], results['maxfreq'])) + log.info(f"\nNumerical dispersion analysis: estimated largest physical \ + phase-velocity error is {results['deltavp']:.2f}% in material \ + '{results['material'].ID}' whose wavelength sampled by {results['N']} \ + cells. Maximum significant frequency estimated as {results['maxfreq']:g}Hz") def reuse_geometry(self): G = self.G # Reset iteration number G.iteration = 0 - self.model_config.inputfilestr = f'\n--- Model {self.model_config.appendmodelnumber}/{self.sim_config.model_end}, input file (not re-processed, i.e. geometry fixed): {self.sim_config.input_file_path}' - log.info(Fore.GREEN + f'{self.model_config.inputfilestr} {'-' * (get_terminal_width() - 1 - len(self.model_config.inputfilestr))}' + Style.RESET_ALL) + self.model_config.inputfilestr = f'\n--- Model {self.model_config.appendmodelnumber}/{self.sim_config.model_end}, \ + input file (not re-processed, i.e. geometry fixed): {self.sim_config.input_file_path}' + log.info(Fore.GREEN + f"{self.model_config.inputfilestr} {'-' * (get_terminal_width() - 1 - len(self.model_config.inputfilestr))}" + Style.RESET_ALL) for grid in [G] + G.subgrids: grid.reset_fields() diff --git a/gprMax/pml.py b/gprMax/pml.py index 7370a759..3fc63ab7 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.py @@ -19,12 +19,11 @@ from importlib import import_module import numpy as np -from tqdm import tqdm import gprMax.config as config -class CFSParameter(object): +class CFSParameter: """Individual CFS parameter (e.g. alpha, kappa, or sigma).""" # Allowable scaling profiles and directions @@ -52,7 +51,7 @@ class CFSParameter(object): self.max = max -class CFS(object): +class CFS: """CFS term for PML.""" def __init__(self): @@ -158,7 +157,7 @@ class CFS(object): return Evalues, Hvalues -class PML(object): +class PML: """Perfectly Matched Layer (PML) Absorbing Boundary Conditions (ABC)""" # Available PML formulations: @@ -177,7 +176,7 @@ class PML(object): def __init__(self, G, ID=None, direction=None, xs=0, xf=0, ys=0, yf=0, zs=0, zf=0): """ Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. ID (str): Identifier for PML slab. direction (str): Direction of increasing absorption. xs, xf, ys, yf, zs, zf (float): Extent of the PML slab. @@ -248,7 +247,7 @@ class PML(object): Args: er (float): Average permittivity of underlying material mr (float): Average permeability of underlying material - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ self.ERA = np.zeros((len(self.CFS), self.thickness), @@ -311,7 +310,7 @@ class PML(object): """This functions updates electric field components with the PML correction. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ pmlmodule = 'gprMax.cython.pml_updates_electric_' + G.pmlformulation @@ -325,7 +324,7 @@ class PML(object): """This functions updates magnetic field components with the PML correction. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ pmlmodule = 'gprMax.cython.pml_updates_magnetic_' + G.pmlformulation @@ -335,16 +334,13 @@ class PML(object): 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. - Args: - G (Grid): Holds essential parameters describing the model. - """ +class CUDAPML(PML): + """Perfectly Matched Layer (PML) Absorbing Boundary Conditions (ABC) for + solving on GPU using CUDA. + """ - self.bpg = (int(np.ceil(((self.EPhi1_gpu.shape[1] + 1) * (self.EPhi1_gpu.shape[2] + 1) * (self.EPhi1_gpu.shape[3] + 1)) / G.tpb[0])), 1, 1) - - def gpu_initialise_arrays(self): + def initialise_arrays(self): """Initialise PML field and coefficient arrays on GPU.""" import pycuda.gpuarray as gpuarray @@ -374,7 +370,16 @@ class PML(object): self.HPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz), dtype=floattype)) self.HPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz), dtype=floattype)) - def gpu_get_update_funcs(self, kernelselectric, kernelsmagnetic): + def set_blocks_per_grid(self, G): + """Set the blocks per grid size used for updating the PML field arrays on a GPU. + + Args: + G (FDTDGrid): Holds essential parameters describing the model. + """ + + self.bpg = (int(np.ceil(((self.EPhi1_gpu.shape[1] + 1) * (self.EPhi1_gpu.shape[2] + 1) * (self.EPhi1_gpu.shape[3] + 1)) / G.tpb[0])), 1, 1) + + def get_update_funcs(self, kernelselectric, kernelsmagnetic): """Get update functions from PML kernels. Args: @@ -387,28 +392,32 @@ class PML(object): self.update_electric_gpu = kernelselectric.get_function('order' + str(len(self.CFS)) + '_' + self.direction) self.update_magnetic_gpu = kernelsmagnetic.get_function('order' + str(len(self.CFS)) + '_' + self.direction) - def gpu_update_electric(self, G): + def update_electric(self, G): """This functions updates electric field components with the PML correction on the GPU. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ self.update_electric_gpu(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_gpu.shape[1]), np.int32(self.EPhi1_gpu.shape[2]), np.int32(self.EPhi1_gpu.shape[3]), np.int32(self.EPhi2_gpu.shape[1]), np.int32(self.EPhi2_gpu.shape[2]), np.int32(self.EPhi2_gpu.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, floattype(self.d), block=G.tpb, grid=self.bpg) - def gpu_update_magnetic(self, G): + def update_magnetic(self, G): """This functions updates magnetic field components with the PML correction on the GPU. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ self.update_magnetic_gpu(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_gpu.shape[1]), np.int32(self.HPhi1_gpu.shape[2]), np.int32(self.HPhi1_gpu.shape[3]), np.int32(self.HPhi2_gpu.shape[1]), np.int32(self.HPhi2_gpu.shape[2]), np.int32(self.HPhi2_gpu.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, floattype(self.d), block=G.tpb, grid=self.bpg) - def pml_information(G): - # No pml + """Information about PMLs. + + Args: + G (FDTDGrid): Holds essential parameters describing the model. + """ + # No PML if all(value == 0 for value in G.pmlthickness.values()): return 'PML: switched off' @@ -417,9 +426,10 @@ def pml_information(G): else: pmlinfo = '' for key, value in G.pmlthickness.items(): - pmlinfo += '{}: {}, '.format(key, value) + pmlinfo += f'{key}: {value}, ' pmlinfo = pmlinfo[:-2] + ' cells' - return '\nPML: formulation: {}, order: {}, thickness: {}'.format(G.pmlformulation, len(G.cfs), pmlinfo) + + return f'\nPML: formulation: {G.pmlformulation}, order: {len(G.cfs)}, thickness: {pmlinfo}' def build_pml(G, key, value): @@ -428,7 +438,7 @@ def build_pml(G, key, value): (based on underlying material er and mr from solid array). Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. key (str): Identifier of PML slab. value (int): Thickness of PML slab in cells. """ diff --git a/gprMax/receivers.py b/gprMax/receivers.py index c7ddda43..4d65b95e 100644 --- a/gprMax/receivers.py +++ b/gprMax/receivers.py @@ -23,7 +23,7 @@ import numpy as np import gprMax.config as config -class Rx(object): +class Rx: """Receiver output points.""" allowableoutputs = ['Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', 'Ix', 'Iy', 'Iz'] @@ -48,7 +48,7 @@ def gpu_initialise_rx_arrays(G): components for receivers. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ import pycuda.gpuarray as gpuarray @@ -79,7 +79,7 @@ def gpu_get_rx_array(rxs_gpu, rxcoords_gpu, G): rxs_gpu (float): numpy array of receiver data from GPU - rows are field components; columns are iterations; pages are receivers. rxcoords_gpu (float): numpy array of receiver coordinates from GPU. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ for rx in G.rxs: diff --git a/gprMax/scene.py b/gprMax/scene.py index 3e99d3a7..44e86108 100644 --- a/gprMax/scene.py +++ b/gprMax/scene.py @@ -34,7 +34,6 @@ class Scene: """Scene stores all of the user created objects.""" def __init__(self): - """Constructor""" self.multiple_cmds = [] self.single_cmds = [] self.geometry_cmds = [] @@ -53,7 +52,7 @@ class Scene: elif isinstance(user_object, UserObjectSingle): self.single_cmds.append(user_object) else: - raise Exception('This Object is Unknown to gprMax') + raise GeneralError('This object is unknown to gprMax') def process_subgrid_commands(self, subgrids): @@ -94,12 +93,12 @@ class Scene: def process_singlecmds(self, G): - # check for duplicate commands and warn user if they exist + # Check for duplicate commands and warn user if they exist cmds_unique = list(set(self.single_cmds)) if len(cmds_unique) != len(self.single_cmds): - raise CmdInputError('Duplicate Single Commands exist in the input.') + raise CmdInputError('Duplicate single-use commands exist in the input.') - # check essential cmds and warn user if missing + # Check essential commands and warn user if missing for cmd_type in self.essential_cmds: d = any([isinstance(cmd, cmd_type) for cmd in cmds_unique]) if not d: diff --git a/gprMax/snapshots.py b/gprMax/snapshots.py index 1caf678d..f27677e0 100644 --- a/gprMax/snapshots.py +++ b/gprMax/snapshots.py @@ -16,9 +16,8 @@ # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -import os -import sys from struct import pack +import sys import numpy as np @@ -27,7 +26,7 @@ from .cython.snapshots import calculate_snapshot_fields from .utilities import round_value -class Snapshot(object): +class Snapshot: """Snapshots of the electric and magnetic field values.""" # Dimensions of largest requested snapshot @@ -151,32 +150,20 @@ class Snapshot(object): self.filehandle = open(self.filename, 'wb') self.filehandle.write('\n'.encode('utf-8')) - self.filehandle.write('\n' - .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')) - 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')) + self.filehandle.write(f'\n'.encode('utf-8')) + self.filehandle.write(f'\n'.encode('utf-8')) + self.filehandle.write(f'\n'.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')) - self.filehandle.write('\n' - .format(Snapshot.floatname, hfield_offset).encode('utf-8')) + self.filehandle.write(f'\n'.encode('utf-8')) + self.filehandle.write(f'\n'.encode('utf-8')) elif self.fieldoutputs['electric']: self.filehandle.write('\n'.encode('utf-8')) - self.filehandle.write('\n' - .format(Snapshot.floatname).encode('utf-8')) + self.filehandle.write(f'\n'.encode('utf-8')) elif self.fieldoutputs['magnetic']: self.filehandle.write('\n'.encode('utf-8')) - self.filehandle.write('\n' - .format(Snapshot.floatname).encode('utf-8')) + self.filehandle.write(f'\n'.encode('utf-8')) self.filehandle.write('\n\n\n\n_'.encode('utf-8')) diff --git a/gprMax/solvers.py b/gprMax/solvers.py index 4e58d7cf..fed2038a 100644 --- a/gprMax/solvers.py +++ b/gprMax/solvers.py @@ -31,7 +31,7 @@ def create_G(sim_config): sim_config (SimConfig): Simulation level configuration object. Returns: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ if sim_config.general['cuda']: @@ -48,7 +48,7 @@ def create_solver(G, sim_config): """Create configured solver object. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. sim_config (SimConfig): simulation level configuration object. Returns: diff --git a/gprMax/sources.py b/gprMax/sources.py index 2067495d..7c2c13f3 100644 --- a/gprMax/sources.py +++ b/gprMax/sources.py @@ -27,7 +27,7 @@ from .grid import Iz from .utilities import round_value -class Source(object): +class Source: """Super-class which describes a generic source.""" def __init__(self): @@ -47,30 +47,33 @@ class Source(object): """Calculates all waveform values for source for duration of simulation. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ # Waveform values for electric sources - calculated half a timestep later - self.waveformvaluesJ = np.zeros((G.iterations), dtype=config.dtypes['float_or_double']) + self.waveformvaluesJ = np.zeros((G.iterations), + dtype=config.dtypes['float_or_double']) # Waveform values for magnetic sources - self.waveformvaluesM = np.zeros((G.iterations), dtype=config.dtypes['float_or_double']) + self.waveformvaluesM = np.zeros((G.iterations), + dtype=config.dtypes['float_or_double']) waveform = next(x for x in G.waveforms if x.ID == self.waveformID) for iteration in range(G.iterations): time = G.dt * iteration if time >= self.start and time <= self.stop: - # Set the time of the waveform evaluation to account for any delay in the start + # Set the time of the waveform evaluation to account for any + # delay in the start time -= self.start self.waveformvaluesJ[iteration] = waveform.calculate_value(time + 0.5 * G.dt, G.dt) self.waveformvaluesM[iteration] = waveform.calculate_value(time, G.dt) class VoltageSource(Source): - """ - A voltage source can be a hard source if it's resistance is zero, i.e. the - time variation of the specified electric field component is prescribed. - If it's resistance is non-zero it behaves as a resistive voltage source. + """A voltage source can be a hard source if it's resistance is zero, + i.e. the time variation of the specified electric field component + is prescribed. If it's resistance is non-zero it behaves as a resistive + voltage source. """ def __init__(self): @@ -82,10 +85,12 @@ class VoltageSource(Source): Args: iteration (int): Current iteration (timestep). - updatecoeffsE (memory view): numpy array of electric field update coefficients. - ID (memory view): numpy array of numeric IDs corresponding to materials in the model. + updatecoeffsE (memory view): numpy array of electric field update + coefficients. + ID (memory view): numpy array of numeric IDs corresponding to + materials in the model. Ex, Ey, Ez (memory view): numpy array of electric field values. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: @@ -96,28 +101,34 @@ class VoltageSource(Source): if self.polarisation == 'x': if self.resistance != 0: - Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * (1 / (self.resistance * G.dy * G.dz)) + Ex[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * + self.waveformvaluesJ[iteration] * + (1 / (self.resistance * G.dy * G.dz))) else: Ex[i, j, k] = -1 * self.waveformvaluesJ[iteration] / G.dx elif self.polarisation == 'y': if self.resistance != 0: - Ey[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * (1 / (self.resistance * G.dx * G.dz)) + Ey[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * + self.waveformvaluesJ[iteration] * + (1 / (self.resistance * G.dx * G.dz))) else: Ey[i, j, k] = -1 * self.waveformvaluesJ[iteration] / G.dy elif self.polarisation == 'z': if self.resistance != 0: - Ez[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * (1 / (self.resistance * G.dx * G.dy)) + Ez[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * + self.waveformvaluesJ[iteration] * + (1 / (self.resistance * G.dx * G.dy))) else: Ez[i, j, k] = -1 * self.waveformvaluesJ[iteration] / G.dz def create_material(self, G): """Create a new material at the voltage source location that adds the - voltage source conductivity to the underlying parameters. + voltage source conductivity to the underlying parameters. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ if self.resistance != 0: @@ -158,10 +169,12 @@ class HertzianDipole(Source): Args: iteration (int): Current iteration (timestep). - updatecoeffsE (memory view): numpy array of electric field update coefficients. - ID (memory view): numpy array of numeric IDs corresponding to materials in the model. + updatecoeffsE (memory view): numpy array of electric field update + coefficients. + ID (memory view): numpy array of numeric IDs corresponding to + materials in the model. Ex, Ey, Ez (memory view): numpy array of electric field values. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: @@ -170,30 +183,35 @@ class HertzianDipole(Source): k = self.zcoord componentID = 'E' + self.polarisation if self.polarisation == 'x': - Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * self.dl * (1 / (G.dx * G.dy * G.dz)) + Ex[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * + self.waveformvaluesJ[iteration] * self.dl * + (1 / (G.dx * G.dy * G.dz))) elif self.polarisation == 'y': - Ey[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * self.dl * (1 / (G.dx * G.dy * G.dz)) + Ey[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * + self.waveformvaluesJ[iteration] * self.dl * + (1 / (G.dx * G.dy * G.dz))) elif self.polarisation == 'z': - Ez[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * self.dl * (1 / (G.dx * G.dy * G.dz)) + Ez[i, j, k] -= (updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * + self.waveformvaluesJ[iteration] * self.dl * + (1 / (G.dx * G.dy * G.dz))) class MagneticDipole(Source): """A magnetic dipole is an additive source (magnetic current density).""" - def __init__(self): - super().__init__() - def update_magnetic(self, iteration, updatecoeffsH, ID, Hx, Hy, Hz, G): """Updates magnetic field values for a magnetic dipole. Args: iteration (int): Current iteration (timestep). - updatecoeffsH (memory view): numpy array of magnetic field update coefficients. - ID (memory view): numpy array of numeric IDs corresponding to materials in the model. + updatecoeffsH (memory view): numpy array of magnetic field update + coefficients. + ID (memory view): numpy array of numeric IDs corresponding to + materials in the model. Hx, Hy, Hz (memory view): numpy array of magnetic field values. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: @@ -203,25 +221,34 @@ class MagneticDipole(Source): componentID = 'H' + self.polarisation if self.polarisation == 'x': - Hx[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesM[iteration] * (1 / (G.dx * G.dy * G.dz)) + Hx[i, j, k] -= (updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * + self.waveformvaluesM[iteration] * + (1 / (G.dx * G.dy * G.dz))) elif self.polarisation == 'y': - Hy[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesM[iteration] * (1 / (G.dx * G.dy * G.dz)) + Hy[i, j, k] -= (updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * + self.waveformvaluesM[iteration] * + (1 / (G.dx * G.dy * G.dz))) elif self.polarisation == 'z': - Hz[i, j, k] -= updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesM[iteration] * (1 / (G.dx * G.dy * G.dz)) + Hz[i, j, k] -= (updatecoeffsH[ID[G.IDlookup[componentID], i, j, k], 4] * + self.waveformvaluesM[iteration] * + (1 / (G.dx * G.dy * G.dz))) def gpu_initialise_src_arrays(sources, G): - """Initialise arrays on GPU for source coordinates/polarisation, other source information, and source waveform values. + """Initialise arrays on GPU for source coordinates/polarisation, other + source information, and source waveform values. Args: sources (list): List of sources of one class, e.g. HertzianDipoles. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. Returns: - srcinfo1_gpu (int): numpy array of source cell coordinates and polarisation information. - srcinfo2_gpu (float): numpy array of other source information, e.g. length, resistance etc... + srcinfo1_gpu (int): numpy array of source cell coordinates and + polarisation information. + srcinfo2_gpu (float): numpy array of other source information, + e.g. length, resistance etc... srcwaves_gpu (float): numpy array of source waveform values. """ @@ -266,7 +293,7 @@ class TransmissionLine(Source): def __init__(self, G): """ Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ super().__init__() @@ -303,7 +330,7 @@ class TransmissionLine(Source): from: http://dx.doi.org/10.1002/mop.10415 Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ for iteration in range(G.iterations): @@ -319,7 +346,7 @@ class TransmissionLine(Source): """Updates absorbing boundary condition at end of the transmission line. Args: - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ h = (config.c * G.dt - self.dl) / (config.c * G.dt + self.dl) @@ -333,11 +360,12 @@ class TransmissionLine(Source): Args: iteration (int): Current iteration (timestep). - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ # Update all the voltage values along the line - self.voltage[1:self.nl] -= self.resistance * (config.c * G.dt / self.dl) * (self.current[1:self.nl] - self.current[0:self.nl - 1]) + self.voltage[1:self.nl] -= (self.resistance * (config.c * G.dt / self.dl) * + (self.current[1:self.nl] - self.current[0:self.nl - 1])) # Update the voltage at the position of the one-way injector excitation self.voltage[self.srcpos] += (config.c * G.dt / self.dl) * self.waveformvaluesJ[iteration] @@ -350,24 +378,30 @@ class TransmissionLine(Source): Args: iteration (int): Current iteration (timestep). - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ # Update all the current values along the line - self.current[0:self.nl - 1] -= (1 / self.resistance) * (config.c * G.dt / self.dl) * (self.voltage[1:self.nl] - self.voltage[0:self.nl - 1]) + self.current[0:self.nl - 1] -= ((1 / self.resistance) * (config.c * G.dt / self.dl) * + (self.voltage[1:self.nl] - self.voltage[0:self.nl - 1])) # Update the current one cell before the position of the one-way injector excitation - self.current[self.srcpos - 1] += (1 / self.resistance) * (config.c * G.dt / self.dl) * self.waveformvaluesM[iteration] + self.current[self.srcpos - 1] += ((1 / self.resistance) * + (config.c * G.dt / self.dl) * + self.waveformvaluesM[iteration]) def update_electric(self, iteration, updatecoeffsE, ID, Ex, Ey, Ez, G): - """Updates electric field value in the main grid from voltage value in the transmission line. + """Updates electric field value in the main grid from voltage value in + the transmission line. Args: iteration (int): Current iteration (timestep). - updatecoeffsE (memory view): numpy array of electric field update coefficients. - ID (memory view): numpy array of numeric IDs corresponding to materials in the model. + updatecoeffsE (memory view): numpy array of electric field update + coefficients. + ID (memory view): numpy array of numeric IDs corresponding to + materials in the model. Ex, Ey, Ez (memory view): numpy array of electric field values. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: @@ -387,14 +421,17 @@ class TransmissionLine(Source): Ez[i, j, k] = - self.voltage[self.antpos] / G.dz def update_magnetic(self, iteration, updatecoeffsH, ID, Hx, Hy, Hz, G): - """Updates current value in transmission line from magnetic field values in the main grid. + """Updates current value in transmission line from magnetic field values + in the main grid. Args: iteration (int): Current iteration (timestep). - updatecoeffsH (memory view): numpy array of magnetic field update coefficients. - ID (memory view): numpy array of numeric IDs corresponding to materials in the model. + updatecoeffsH (memory view): numpy array of magnetic field update + coefficients. + ID (memory view): numpy array of numeric IDs corresponding to + materials in the model. Hx, Hy, Hz (memory view): numpy array of magnetic field values. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing the model. """ if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: diff --git a/gprMax/subgrids/updates.py b/gprMax/subgrids/updates.py index 77c32e3c..ae77c0d4 100644 --- a/gprMax/subgrids/updates.py +++ b/gprMax/subgrids/updates.py @@ -75,7 +75,7 @@ class SubgridUpdater(CPUUpdates): subgrid (SubGrid3d): Subgrid to be updated. precursors (PrecursorNodes): Precursor nodes associated with the subgrid - contain interpolated fields. - G (Grid): Holds essential parameters describing the model. + G (FDTDGrid): Holds essential parameters describing a model. """ super().__init__(subgrid) self.precursors = precursors diff --git a/gprMax/updates.py b/gprMax/updates.py index a30ac8a0..489ef16a 100644 --- a/gprMax/updates.py +++ b/gprMax/updates.py @@ -35,16 +35,16 @@ from .snapshots import gpu_get_snapshot_array from .sources import gpu_initialise_src_arrays from .utilities import timer - log = logging.getLogger(__name__) + class CPUUpdates: """Defines update functions for CPU-based solver.""" def __init__(self, G): """ Args: - G (FDTDGrid): FDTD grid object + G (FDTDGrid): Holds essential parameters describing the model. """ self.grid = G self.dispersive_update_a = None @@ -143,14 +143,21 @@ class CPUUpdates: update any Hertzian dipole sources last. """ for source in self.grid.voltagesources + self.grid.transmissionlines + self.grid.hertziandipoles: - source.update_electric(self.grid.iteration, self.grid.updatecoeffsE, self.grid.ID, self.grid.Ex, self.grid.Ey, self.grid.Ez, self.grid) + source.update_electric(self.grid.iteration, + self.grid.updatecoeffsE, + self.grid.ID, + self.grid.Ex, + self.grid.Ey, + self.grid.Ez, + self.grid) self.grid.iteration += 1 def update_electric_b(self): - """If there are any dispersive materials do 2nd part of dispersive update - - it is split into two parts as it requires present and updated electric - field values. Therefore it can only be completely updated after the - electric field has been updated by the PML and source updates. + """If there are any dispersive materials do 2nd part of dispersive + update - it is split into two parts as it requires present and + updated electric field values. Therefore it can only be completely + updated after the electric field has been updated by the PML and + source updates. """ if config.materials['maxpoles'] != 0: self.dispersive_update_b(self.grid.nx, @@ -236,7 +243,7 @@ class CUDAUpdates: def __init__(self, G): """ Args: - G (FDTDGrid): FDTD grid object + G (FDTDGrid): Holds essential parameters describing the model. """ self.grid = G @@ -248,6 +255,7 @@ class CUDAUpdates: drv.init() # Suppress nvcc warnings on Windows + log.debug('Move nvcc compiler options to simulation config') if sys.platform == 'win32': self.compiler_opts = ['-w'] else: @@ -269,9 +277,41 @@ class CUDAUpdates: get kernel functions. """ if config.materials['maxpoles'] > 0: - kernels_fields = SourceModule(kernels_template_fields.substitute(REAL=cudafloattype, COMPLEX=cudacomplextype, N_updatecoeffsE=self.grid.updatecoeffsE.size, N_updatecoeffsH=self.grid.updatecoeffsH.size, NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], NY_MATDISPCOEFFS=self.grid.updatecoeffsdispersive.shape[1], NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1, NX_ID=self.grid.ID.shape[1], NY_ID=self.grid.ID.shape[2], NZ_ID=self.grid.ID.shape[3], NX_T=self.grid.Tx.shape[1], NY_T=self.grid.Tx.shape[2], NZ_T=self.grid.Tx.shape[3]), options=self.compiler_opts) + kernels_fields = SourceModule(kernels_template_fields.substitute( + REAL=cudafloattype, + COMPLEX=cudacomplextype, + N_updatecoeffsE=self.grid.updatecoeffsE.size, + N_updatecoeffsH=self.grid.updatecoeffsH.size, + NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], + NY_MATDISPCOEFFS=self.grid.updatecoeffsdispersive.shape[1], + NX_FIELDS=self.grid.nx + 1, + NY_FIELDS=self.grid.ny + 1, + NZ_FIELDS=self.grid.nz + 1, + NX_ID=self.grid.ID.shape[1], + NY_ID=self.grid.ID.shape[2], + NZ_ID=self.grid.ID.shape[3], + NX_T=self.grid.Tx.shape[1], + NY_T=self.grid.Tx.shape[2], + NZ_T=self.grid.Tx.shape[3]), + options=self.compiler_opts) else: # Set to one any substitutions for dispersive materials - kernels_fields = SourceModule(kernels_template_fields.substitute(REAL=cudafloattype, COMPLEX=cudacomplextype, N_updatecoeffsE=self.grid.updatecoeffsE.size, N_updatecoeffsH=self.grid.updatecoeffsH.size, NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], NY_MATDISPCOEFFS=1, NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1, NX_ID=self.grid.ID.shape[1], NY_ID=self.grid.ID.shape[2], NZ_ID=self.grid.ID.shape[3], NX_T=1, NY_T=1, NZ_T=1), options=self.compiler_opts) + kernels_fields = SourceModule(kernels_template_fields.substitute( + REAL=cudafloattype, + COMPLEX=cudacomplextype, + N_updatecoeffsE=self.grid.updatecoeffsE.size, + N_updatecoeffsH=self.grid.updatecoeffsH.size, + NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], + NY_MATDISPCOEFFS=1, + NX_FIELDS=self.grid.nx + 1, + NY_FIELDS=self.grid.ny + 1, + NZ_FIELDS=self.grid.nz + 1, + NX_ID=self.grid.ID.shape[1], + NY_ID=self.grid.ID.shape[2], + NZ_ID=self.grid.ID.shape[3], + NX_T=1, + NY_T=1, + NZ_T=1), + options=self.compiler_opts) self.update_electric = kernels_fields.get_function("update_electric") self.update_magnetic = kernels_fields.get_function("update_magnetic") if self.grid.updatecoeffsE.nbytes + self.grid.updatecoeffsH.nbytes > self.grid.gpu.constmem: @@ -292,11 +332,35 @@ class CUDAUpdates: """PMLS - prepare kernels and get kernel functions.""" if self.grid.pmls: pmlmodulelectric = 'gprMax.cuda.pml_updates_electric_' + self.grid.pmlformulation - kernelelectricfunc = getattr(import_module(pmlmodulelectric), 'kernels_template_pml_electric_' + self.grid.pmlformulation) + kernelelectricfunc = getattr(import_module(pmlmodulelectric), + 'kernels_template_pml_electric_' + + self.grid.pmlformulation) pmlmodulemagnetic = 'gprMax.cuda.pml_updates_magnetic_' + self.grid.pmlformulation - kernelmagneticfunc = getattr(import_module(pmlmodulemagnetic), 'kernels_template_pml_magnetic_' + self.grid.pmlformulation) - kernels_pml_electric = SourceModule(kernelelectricfunc.substitute(REAL=cudafloattype, N_updatecoeffsE=self.grid.updatecoeffsE.size, NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1, NX_ID=self.grid.ID.shape[1], NY_ID=self.grid.ID.shape[2], NZ_ID=self.grid.ID.shape[3]), options=self.compiler_opts) - kernels_pml_magnetic = SourceModule(kernelmagneticfunc.substitute(REAL=cudafloattype, N_updatecoeffsH=self.grid.updatecoeffsH.size, NY_MATCOEFFS=self.grid.updatecoeffsH.shape[1], NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1, NX_ID=self.gridG.ID.shape[1], NY_ID=self.grid.ID.shape[2], NZ_ID=self.grid.ID.shape[3]), options=self.compiler_opts) + kernelmagneticfunc = getattr(import_module(pmlmodulemagnetic), + 'kernels_template_pml_magnetic_' + + self.grid.pmlformulation) + kernels_pml_electric = SourceModule(kernelelectricfunc.substitute( + REAL=cudafloattype, + N_updatecoeffsE=self.grid.updatecoeffsE.size, + NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], + NX_FIELDS=self.grid.nx + 1, + NY_FIELDS=self.grid.ny + 1, + NZ_FIELDS=self.grid.nz + 1, + NX_ID=self.grid.ID.shape[1], + NY_ID=self.grid.ID.shape[2], + NZ_ID=self.grid.ID.shape[3]), + options=self.compiler_opts) + kernels_pml_magnetic = SourceModule(kernelmagneticfunc.substitute( + REAL=cudafloattype, + N_updatecoeffsH=self.grid.updatecoeffsH.size, + NY_MATCOEFFS=self.grid.updatecoeffsH.shape[1], + NX_FIELDS=self.grid.nx + 1, + NY_FIELDS=self.grid.ny + 1, + NZ_FIELDS=self.grid.nz + 1, + NX_ID=self.gridG.ID.shape[1], + NY_ID=self.grid.ID.shape[2], + NZ_ID=self.grid.ID.shape[3]), + options=self.compiler_opts) self.copy_mat_coeffs() # Set block per grid, initialise arrays on GPU, and get kernel functions for pml in self.grid.pmls: @@ -305,20 +369,42 @@ class CUDAUpdates: pml.gpu_set_blocks_per_grid(self.grid) def set_rx_kernel(self): - """Receivers - initialise arrays on GPU, prepare kernel and - get kernel function. + """Receivers - initialise arrays on GPU, prepare kernel and get kernel + function. """ if self.grid.rxs: rxcoords_gpu, rxs_gpu = gpu_initialise_rx_arrays(self.grid) - kernel_store_outputs = SourceModule(kernel_template_store_outputs.substitute(REAL=cudafloattype, NY_RXCOORDS=3, NX_RXS=6, NY_RXS=self.grid.iterations, NZ_RXS=len(self.grid.rxs), NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1), options=self.compiler_opts) + kernel_store_outputs = SourceModule(kernel_template_store_outputs.substitute( + REAL=cudafloattype, + NY_RXCOORDS=3, + NX_RXS=6, + NY_RXS=self.grid.iterations, + NZ_RXS=len(self.grid.rxs), + NX_FIELDS=self.grid.nx + 1, + NY_FIELDS=self.grid.ny + 1, + NZ_FIELDS=self.grid.nz + 1), + options=self.compiler_opts) self.store_outputs = kernel_store_outputs.get_function("store_outputs") def set_src_kernels(self): - """Sources - initialise arrays on GPU, prepare kernel and - get kernel function. + """Sources - initialise arrays on GPU, prepare kernel and get kernel + function. """ if self.grid.voltagesources + self.grid.hertziandipoles + self.grid.magneticdipoles: - kernels_sources = SourceModule(kernels_template_sources.substitute(REAL=cudafloattype, N_updatecoeffsE=self.grid.updatecoeffsE.size, N_updatecoeffsH=self.grid.updatecoeffsH.size, NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], NY_SRCINFO=4, NY_SRCWAVES=self.grid.iterations, NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1, NX_ID=self.grid.ID.shape[1], NY_ID=self.grid.ID.shape[2], NZ_ID=self.grid.ID.shape[3]), options=self.compiler_opts) + kernels_sources = SourceModule(kernels_template_sources.substitute( + REAL=cudafloattype, + N_updatecoeffsE=self.grid.updatecoeffsE.size, + N_updatecoeffsH=self.grid.updatecoeffsH.size, + NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], + NY_SRCINFO=4, + NY_SRCWAVES=self.grid.iterations, + NX_FIELDS=self.grid.nx + 1, + NY_FIELDS=self.grid.ny + 1, + NZ_FIELDS=self.grid.nz + 1, + NX_ID=self.grid.ID.shape[1], + NY_ID=self.grid.ID.shape[2], + NZ_ID=self.grid.ID.shape[3]), + options=self.compiler_opts) self.copy_mat_coeffs() if self.grid.hertziandipoles: self.srcinfo1_hertzian_gpu, self.srcinfo2_hertzian_gpu, self.srcwaves_hertzian_gpu = gpu_initialise_src_arrays(self.grid.hertziandipoles, self.grid) @@ -331,12 +417,20 @@ class CUDAUpdates: self.update_voltage_source_gpu = kernels_sources.get_function("update_voltage_source") def set_snapshot_kernel(self): - """Snapshots - initialise arrays on GPU, prepare kernel and - get kernel function. + """Snapshots - initialise arrays on GPU, prepare kernel and get kernel + function. """ if self.grid.snapshots: self.snapEx_gpu, self.snapEy_gpu, self.snapEz_gpu, self.snapHx_gpu, self.snapHy_gpu, self.snapHz_gpu = gpu_initialise_snapshot_array(self.grid) - kernel_store_snapshot = SourceModule(kernel_template_store_snapshot.substitute(REAL=cudafloattype, NX_SNAPS=Snapshot.nx_max, NY_SNAPS=Snapshot.ny_max, NZ_SNAPS=Snapshot.nz_max, NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1), options=self.compiler_opts) + kernel_store_snapshot = SourceModule(kernel_template_store_snapshot.substitute( + REAL=cudafloattype, + NX_SNAPS=Snapshot.nx_max, + NY_SNAPS=Snapshot.ny_max, + NZ_SNAPS=Snapshot.nz_max, + NX_FIELDS=self.grid.nx + 1, + NY_FIELDS=self.grid.ny + 1, + NZ_FIELDS=self.grid.nz + 1), + options=self.compiler_opts) self.store_snapshot_gpu = kernel_store_snapshot.get_function("store_snapshot") def copy_mat_coeffs(self): @@ -528,10 +622,11 @@ class CUDAUpdates: self.grid.iteration += 1 def update_electric_b(self): - """If there are any dispersive materials do 2nd part of dispersive update - - it is split into two parts as it requires present and updated electric - field values. Therefore it can only be completely updated after the - electric field has been updated by the PML and source updates. + """If there are any dispersive materials do 2nd part of dispersive + update - it is split into two parts as it requires present and + updated electric field values. Therefore it can only be completely + updated after the electric field has been updated by the PML and + source updates. """ if config.materials['maxpoles'] != 0: self.dispersive_update_b(np.int32(self.grid.nx), diff --git a/gprMax/user_inputs.py b/gprMax/user_inputs.py index 67a93fd9..85b624bd 100644 --- a/gprMax/user_inputs.py +++ b/gprMax/user_inputs.py @@ -16,6 +16,8 @@ # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . +import logging + from colorama import init from colorama import Fore from colorama import Style @@ -27,20 +29,22 @@ from .exceptions import CmdInputError from .subgrids.base import SubGridBase from .utilities import round_value +log = logging.getLogger(__name__) """Module contains classes to handle points supplied by a user. The -classes implement a common interface such that geometry building objects -such as box or triangle do not need to have any knowledge which grid to which -they are rounding continuous points or checking the point is within the grid. -Additionally all logic related to rounding points etc is encapulsated here. + classes implement a common interface such that geometry building objects + such as box or triangle do not need to have any knowledge which grid to + which they are rounding continuous points or checking the point is within + the grid. Additionally all logic related to rounding points etc is + encapulsated here. """ def create_user_input_points(grid, user_obj): """Return a point checker class based on the grid supplied.""" - if isinstance(grid, SubGridBase): - # local object configuration trumps. User can turn of autotranslate for + if isinstance(grid, SubGridBase): + # Local object configuration trumps. User can turn of autotranslate for # specfic objects. if not user_obj.autotranslate and config.general['autotranslate']: return MainGridUserInput(grid) @@ -69,12 +73,12 @@ class UserInput: # Incorrect index i = p[v.index(err.args[0])] if name: - fs = "'{}' the {} {}-coordinate {:g} is not within the model domain" - s = fs.format(cmd_str, err.args[0], name, i * dl) + s = f"'{cmd_str}' the {err.args[0]} {name}-coordinate {i * dl:g} \ + is not within the model domain" else: - fs = "'{}' {}-coordinate {:g} is not within the model domain" - s = fs.format(cmd_str, err.args[0], i * dl) - raise CmdInputError(s) + s = f"'{cmd_str}' {err.args[0]}-coordinate {i * dl:g} is not \ + within the model domain" + raise CmdInputError(log.exception(s)) def discretise_point(self, p): """Function to get the index of a continuous point with the grid.""" @@ -102,10 +106,8 @@ class MainGridUserInput(UserInput): p = self.check_point(p, cmd_str, name) if self.grid.within_pml(p): - s = """WARNING: '{}' sources and receivers should not normally be \ - positioned within the PML.""" - fs = s.format(cmd_str) - print(Fore.RED + fs + Style.RESET_ALL) + log.warning(Fore.RED + f"'{cmd_str}' sources and receivers should \ + not normally be positioned within the PML." + Style.RESET_ALL) return p @@ -114,9 +116,8 @@ class MainGridUserInput(UserInput): p2 = self.check_point(p2, cmd_str, name='upper') if np.greater(p1, p2).any(): - s = """'{}' the lower coordinates should be less than the upper \ - coordinates""" - raise CmdInputError(s.format(cmd_str)) + raise CmdInputError(log.exception(f"'{cmd_str}' the lower coordinates \ + should be less than the upper coordinates.")) return p1, p2 @@ -130,14 +131,14 @@ class MainGridUserInput(UserInput): class SubgridUserInput(MainGridUserInput): """Class to handle (x, y, z) points supplied by the user in the sub grid. - This class autotranslates points from main grid to subgrid equivalent (within IS). - Useful if material traverse is not required. + This class autotranslates points from main grid to subgrid equivalent + (within IS). Useful if material traverse is not required. """ def __init__(self, grid): super().__init__(grid) - # defines the region exposed to the user + # Defines the region exposed to the user self.inner_bound = np.array([ grid.n_boundary_cells_x, grid.n_boundary_cells_y, @@ -148,7 +149,8 @@ class SubgridUserInput(MainGridUserInput): def translate_to_gap(self, p): """Function to translate the user input point to the real point in the - subgrid""" + subgrid. + """ p1 = (p[0] - self.grid.i0 * self.grid.ratio) + self.grid.n_boundary_cells_x p2 = (p[1] - self.grid.j0 * self.grid.ratio) + self.grid.n_boundary_cells_y @@ -157,9 +159,10 @@ class SubgridUserInput(MainGridUserInput): return np.array([p1, p2, p3]) def discretise_point(self, p): - """"Function to discretise a point. Does not provide any checks. The - user enters coordinates relative to self.inner_bound. This function - translate the user point to the correct index for building objects""" + """Function to discretise a point. Does not provide any checks. The + user enters coordinates relative to self.inner_bound. This function + translate the user point to the correct index for building objects. + """ p = super().discretise_point(p) p_t = self.translate_to_gap(p) @@ -177,8 +180,6 @@ class SubgridUserInput(MainGridUserInput): # the OS non-working region. if (np.less(p_t, self.inner_bound).any() or np.greater(p_t, self.outer_bound).any()): - s = """WARNING: '{}' This object traverses the Outer Surface. This \ - is an advanced feature.""" - print(Fore.RED + s.format(cmd_str) + Style.RESET_ALL) - + log.warning(Fore.RED + f"'{cmd_str}' this object traverses the \ + Outer Surface. This is an advanced feature." + Style.RESET_ALL) return p_t diff --git a/gprMax/utilities.py b/gprMax/utilities.py index 6f2fdfb5..cc06de60 100644 --- a/gprMax/utilities.py +++ b/gprMax/utilities.py @@ -19,6 +19,7 @@ from contextlib import contextmanager import codecs import decimal as d +import logging import os import platform import psutil @@ -37,6 +38,8 @@ import numpy as np from .exceptions import GeneralError +log = logging.getLogger(__name__) + def get_terminal_width(): """Get/set width of terminal being used. @@ -74,7 +77,7 @@ def logo(version): |___/|_| v""" + version - log.info(f'{description} {'=' * (get_terminal_width() - len(description) - 1)}\n') + log.info(f"{description} {'=' * (get_terminal_width() - len(description) - 1)}\n") log.info(Fore.CYAN + f'{logo}\n') log.info(Style.RESET_ALL + textwrap.fill(copyright, width=get_terminal_width() - 1, initial_indent=' ')) log.info(textwrap.fill(authors, width=get_terminal_width() - 1, initial_indent=' ')) @@ -336,7 +339,44 @@ def get_host_info(): return hostinfo -class GPU(object): +def set_openmp_threads(): + """Sets the number of OpenMP CPU threads for parallelised parts of code. + + Returns: + nthreads (int): Number of OpenMP threads. + """ + + if sys.platform == 'darwin': + # Should waiting threads consume CPU power (can drastically effect + # performance) + os.environ['OMP_WAIT_POLICY'] = 'ACTIVE' + # Number of threads may be adjusted by the run time environment to best + # utilize system resources + os.environ['OMP_DYNAMIC'] = 'FALSE' + # Each place corresponds to a single core (having one or more hardware threads) + os.environ['OMP_PLACES'] = 'cores' + # Bind threads to physical cores + os.environ['OMP_PROC_BIND'] = 'TRUE' + # Prints OMP version and environment variables (useful for debug) + # os.environ['OMP_DISPLAY_ENV'] = 'TRUE' + + # Catch bug with Windows Subsystem for Linux (https://github.com/Microsoft/BashOnWindows/issues/785) + if 'Microsoft' in config.hostinfo['osversion']: + os.environ['KMP_AFFINITY'] = 'disabled' + del os.environ['OMP_PLACES'] + del os.environ['OMP_PROC_BIND'] + + if os.environ.get('OMP_NUM_THREADS'): + nthreads = int(os.environ.get('OMP_NUM_THREADS')) + else: + # Set number of threads to number of physical CPU cores + nthreads = config.hostinfo['physicalcores'] + os.environ['OMP_NUM_THREADS'] = str(G.nthreads) + + return nthreads + + +class GPU: """GPU information.""" def __init__(self, deviceID): diff --git a/gprMax/waveforms.py b/gprMax/waveforms.py index 1da9f262..667b4114 100644 --- a/gprMax/waveforms.py +++ b/gprMax/waveforms.py @@ -16,21 +16,29 @@ # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . +import logging + import numpy as np +log = logging.getLogger(__name__) -class Waveform(object): + +class Waveform: """Definitions of waveform shapes that can be used with sources.""" - types = ['gaussian', 'gaussiandot', 'gaussiandotnorm', 'gaussiandotdot', 'gaussiandotdotnorm', 'gaussianprime', 'gaussiandoubleprime', 'ricker', 'sine', 'contsine', 'impulse', 'user'] + types = ['gaussian', 'gaussiandot', 'gaussiandotnorm', 'gaussiandotdot', + 'gaussiandotdotnorm', 'gaussianprime', 'gaussiandoubleprime', + 'ricker', 'sine', 'contsine', 'impulse', 'user'] # Information about specific waveforms: # - # gaussianprime and gaussiandoubleprime waveforms are the first derivative and second derivative of the 'base' gaussian - # waveform, i.e. the centre frequencies of the waveforms will rise for the first and second derivatives. + # gaussianprime and gaussiandoubleprime waveforms are the first derivative + # and second derivative of the 'base' gaussian waveform, i.e. the centre + # frequencies of the waveforms will rise for the first and second derivatives. # - # gaussiandot, gaussiandotnorm, gaussiandotdot, gaussiandotdotnorm, ricker waveforms have their centre frequencies - # specified by the user, i.e. they are not derived from the 'base' gaussian + # gaussiandot, gaussiandotnorm, gaussiandotdot, gaussiandotdotnorm, + # ricker waveforms have their centre frequencies specified by the user, + # i.e. they are not derived from the 'base' gaussian def __init__(self): self.ID = None @@ -43,7 +51,9 @@ class Waveform(object): self.delay = 0 def calculate_coefficients(self): - """Calculates coefficients (used to calculate values) for specific waveforms.""" + """Calculates coefficients (used to calculate values) for specific + waveforms. + """ if self.type == 'gaussian' or self.type == 'gaussiandot' or self.type == 'gaussiandotnorm' or self.type == 'gaussianprime' or self.type == 'gaussiandoubleprime': self.chi = 1 / self.freq @@ -81,17 +91,20 @@ class Waveform(object): elif self.type == 'gaussiandotdot' or self.type == 'gaussiandoubleprime': delay = time - self.chi - ampvalue = 2 * self.zeta * (2 * self.zeta * delay**2 - 1) * np.exp(-self.zeta * delay**2) + ampvalue = (2 * self.zeta * (2 * self.zeta * delay**2 - 1) * + np.exp(-self.zeta * delay**2)) elif self.type == 'gaussiandotdotnorm': delay = time - self.chi normalise = 1 / (2 * self.zeta) - ampvalue = 2 * self.zeta * (2 * self.zeta * delay**2 - 1) * np.exp(-self.zeta * delay**2) * normalise + ampvalue = (2 * self.zeta * (2 * self.zeta * delay**2 - 1) * + np.exp(-self.zeta * delay**2) * normalise) elif self.type == 'ricker': delay = time - self.chi normalise = 1 / (2 * self.zeta) - ampvalue = - (2 * self.zeta * (2 * self.zeta * delay**2 - 1) * np.exp(-self.zeta * delay**2)) * normalise + ampvalue = - ((2 * self.zeta * (2 * self.zeta * delay**2 - 1) * + np.exp(-self.zeta * delay**2)) * normalise) elif self.type == 'sine': ampvalue = np.sin(2 * np.pi * self.freq * time) @@ -120,6 +133,5 @@ class Waveform(object): return ampvalue def __str__(self): - fs = 'Waveform: ID={}, type={}, amp{}, freq={}' - s = fs.format(self.ID, self.type, self.amp, self.freq) - return s + log.debug('Do we need this?') + return f'Waveform: ID={self.ID}, type={self.type}, amp{self.amp}, freq={self.freq}'