publicImage/gprMax
1
0
派生 0
镜像自地址 https://gitee.com/sunhf/gprMax.git 已同步 2025-08-08 07:24:19 +08:00
代码 工单 软件包 项目 版本发布 百科 动态

More changes discussed with John. F-strings, logging etc...

浏览代码
这个提交包含在：
Craig Warren
2019-10-11 11:08:51 +01:00
父节点 c11dbd7671
当前提交 1ad97fbeef
共有 26 个文件被更改，包括 550 次插入 和 426 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

2
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

108
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

71
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):

20
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

7
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 <http://www.gnu.org/licenses/>.
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)

8
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')

6
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):

15
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('</gprMax>\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'

26
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:

28
gprMax/input_cmds_file.py → 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 <http://www.gnu.org/licenses/>.
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

8
gprMax/input_cmds_geometry.py → 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 <http://www.gnu.org/licenses/>.
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

37
gprMax/input_cmds_multiuse.py → 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 <http://www.gnu.org/licenses/>.
# 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 <http://www.gnu.org/licenses/>.
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

23
gprMax/input_cmds_singleuse.py → 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 <http://www.gnu.org/licenses/>.
# 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 <http://www.gnu.org/licenses/>.
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)

16
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)

56
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()

62
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.
"""

6
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:

9
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:

31
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 <http://www.gnu.org/licenses/>.
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('<?xml version="1.0"?>\n'.encode('utf-8'))
self.filehandle.write('<VTKFile type="ImageData" version="1.0" byte_order="{}">\n'
.format(Snapshot.byteorder).encode('utf-8'))
self.filehandle.write('<ImageData WholeExtent="{} {} {} {} {} {}" Origin="0 0 0" Spacing="{:.3} {:.3} {:.3}">\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('<Piece Extent="{} {} {} {} {} {}">\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'<VTKFile type="ImageData" version="1.0" byte_order="{Snapshot.byteorder}">\n'.encode('utf-8'))
self.filehandle.write(f'<ImageData WholeExtent="{self.xs} {round_value(self.xf / self.dx)} {self.ys} {round_value(self.yf / self.dy)} {self.zs} {round_value(self.zf / self.dz)}" Origin="0 0 0" Spacing="{self.dx * G.dx:.3} {self.dy * G.dy:.3} {self.dz * G.dz:.3}">\n'.encode('utf-8'))
self.filehandle.write(f'<Piece Extent="{self.xs} {round_value(self.xf / self.dx)} {self.ys} {round_value(self.yf / self.dy)} {self.zs} {round_value(self.zf / self.dz)}">\n'.encode('utf-8'))
if self.fieldoutputs['electric'] and self.fieldoutputs['magnetic']:
self.filehandle.write('<CellData Vectors="E-field H-field">\n'.encode('utf-8'))
self.filehandle.write('<DataArray type="{}" Name="E-field" NumberOfComponents="3" format="appended" offset="0" />\n'
.format(Snapshot.floatname).encode('utf-8'))
self.filehandle.write('<DataArray type="{}" Name="H-field" NumberOfComponents="3" format="appended" offset="{}" />\n'
.format(Snapshot.floatname, hfield_offset).encode('utf-8'))
self.filehandle.write(f'<DataArray type="{Snapshot.floatname}" Name="E-field" NumberOfComponents="3" format="appended" offset="0" />\n'.encode('utf-8'))
self.filehandle.write(f'<DataArray type="{Snapshot.floatname}" Name="H-field" NumberOfComponents="3" format="appended" offset="{hfield_offset}" />\n'.encode('utf-8'))
elif self.fieldoutputs['electric']:
self.filehandle.write('<CellData Vectors="E-field">\n'.encode('utf-8'))
self.filehandle.write('<DataArray type="{}" Name="E-field" NumberOfComponents="3" format="appended" offset="0" />\n'
.format(Snapshot.floatname).encode('utf-8'))
self.filehandle.write(f'<DataArray type="{Snapshot.floatname}" Name="E-field" NumberOfComponents="3" format="appended" offset="0" />\n'.encode('utf-8'))
elif self.fieldoutputs['magnetic']:
self.filehandle.write('<CellData Vectors="H-field">\n'.encode('utf-8'))
self.filehandle.write('<DataArray type="{}" Name="H-field" NumberOfComponents="3" format="appended" offset="0" />\n'
.format(Snapshot.floatname).encode('utf-8'))
self.filehandle.write(f'<DataArray type="{Snapshot.floatname}" Name="H-field" NumberOfComponents="3" format="appended" offset="0" />\n'.encode('utf-8'))
self.filehandle.write('</CellData>\n</Piece>\n</ImageData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))

4
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:

141
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:

2
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

149
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),

59
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 <http://www.gnu.org/licenses/>.
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

44
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):

38
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 <http://www.gnu.org/licenses/>.
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}'