From 90abd52ce4b9cc87302efd114dcf2ce5aa6955aa Mon Sep 17 00:00:00 2001 From: John Hartley Date: Mon, 5 Aug 2019 21:50:49 +0100 Subject: [PATCH] working up to solver.solve() --- gprMax/cmds_single_use.py | 29 +++++++++++------------ gprMax/config.py | 3 ++- gprMax/grid.py | 50 ++++++++++++++++++++++++++++++++------- gprMax/model_build_run.py | 6 +++-- gprMax/pml.py | 2 +- gprMax/scene.py | 2 +- gprMax/user_inputs.py | 1 + 7 files changed, 64 insertions(+), 29 deletions(-) diff --git a/gprMax/cmds_single_use.py b/gprMax/cmds_single_use.py index a856c858..d0ebfdc6 100644 --- a/gprMax/cmds_single_use.py +++ b/gprMax/cmds_single_use.py @@ -19,15 +19,17 @@ import inspect import os import sys +import decimal as d from .config import c from .config import dtypes +from .config import hostinfo from .exceptions import CmdInputError from .waveforms import Waveform +from .utilities import round_value import numpy as np from scipy import interpolate - from colorama import init from colorama import Fore from colorama import Style @@ -107,16 +109,14 @@ class DomainSingle(UserObjectSingle): G.mode = '3D' # Round down time step to nearest float with precision one less than hardware maximum. Avoids inadvertently exceeding the CFL due to binary representation of floating point number. - G.round_time_step() + # Round down time step to nearest float with precision one less than hardware maximum. + # Avoids inadvertently exceeding the CFL due to binary representation of floating point number. + G.dt = round_value(G.dt, decimalplaces=d.getcontext().prec - 1) if G.messages: print('Mode: {}'.format(G.mode)) print('Time step (at CFL limit): {:g} secs'.format(G.dt)) - # do threads here - from .utilities import get_host_info - hostinfo = get_host_info() - # 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) @@ -126,7 +126,7 @@ class DomainSingle(UserObjectSingle): # 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 G.hostinfo['osversion']: + if 'Microsoft' in hostinfo['osversion']: os.environ['KMP_AFFINITY'] = 'disabled' del os.environ['OMP_PLACES'] del os.environ['OMP_PROC_BIND'] @@ -140,7 +140,7 @@ class DomainSingle(UserObjectSingle): if G.messages: print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads)) - if G.nthreads > G.hostinfo['physicalcores']: + if G.nthreads > hostinfo['physicalcores']: print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, hostinfo['physicalcores']) + Style.RESET_ALL) @@ -180,20 +180,21 @@ class Discretisation(UserObjectSingle): def create(self, G, uip): try: + G.dl = np.array(self.kwargs['p1']) G.dx, G.dy, G.dz = self.kwargs['p1'] except KeyError: raise CmdInputError('Discretisation requires a point') - if G.dx <= 0: + if G.dl[0] <= 0: raise CmdInputError('Discretisation requires the x-direction spatial step to be greater than zero') - if G.dy <= 0: + if G.dl[1] <= 0: raise CmdInputError(' Discretisation requires the y-direction spatial step to be greater than zero') - if G.dz <= 0: + if G.dl[2] <= 0: raise CmdInputError('Discretisation requires the z-direction spatial step to be greater than zero') if G.messages: - print('Spatial discretisation: {:g} x {:g} x {:g}m'.format(G.dx, G.dy, G.dz)) + print('Spatial discretisation: {:g} x {:g} x {:g}m'.format(*G.dl)) class TimeWindow(UserObjectSingle): @@ -307,8 +308,6 @@ class NumThreads(UserObjectSingle): def create(self, G, uip): # Get information about host machine - from .utilities import get_host_info - hostinfo = get_host_info() try: n = self.kwargs['n'] @@ -323,7 +322,7 @@ class NumThreads(UserObjectSingle): if G.messages: print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads)) - if G.nthreads > G.hostinfo['physicalcores']: + if G.nthreads > hostinfo['physicalcores']: print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, hostinfo['physicalcores']) + Style.RESET_ALL) # Print information about any GPU in use diff --git a/gprMax/config.py b/gprMax/config.py index 5a35c55a..d162da63 100644 --- a/gprMax/config.py +++ b/gprMax/config.py @@ -46,7 +46,7 @@ z0 = np.sqrt(m0 / e0) # mode: 2D TMx, 2D TMy, 2D TMz, or 3D # cpu, cuda, opencl: solver type general = {'inputfilepath': '', 'outputfilepath': '', 'messages': True, - 'progressbars': True, 'mode': None, 'cpu': True, 'cuda': False, 'opencl': False} + 'progressbars': True, 'mode': '3D', 'cpu': True, 'cuda': False, 'opencl': False} # Store information about host machine hostinfo = get_host_info() @@ -141,6 +141,7 @@ class SimulationConfig: self.general = {} self.general['messages'] = True self.geometry_fixed = args.geometry_fixed + self.geometry_only = args.geometry_only self.scenes = args.scenes self.set_input_file_path() diff --git a/gprMax/grid.py b/gprMax/grid.py index 8de05811..411d4272 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -65,14 +65,38 @@ class Grid(object): def get(self, i, j, k): return self.grid[i, j, k] - def within_bounds(self, **kwargs): - for co, val in kwargs.items(): - if val < 0 or val > getattr(self, 'n' + co): - raise ValueError(co) + def within_bounds(self, p): + if p[0] < 0 or p[0] > self.nx: + raise ValueError('x') + if p[1] < 0 or p[1] > self.ny: + raise ValueError('y') + if p[2] < 0 or p[2] > self.nz: + raise ValueError('z') - def calculate_coord(self, coord, val): - co = round_value(float(val) / getattr(self, 'd' + coord)) - return co + def discretise_point(self, p): + x = round_value(float(p[0]) / self.dx) + y = round_value(float(p[1]) / self.dy) + z = round_value(float(p[2]) / self.dz) + return (x, y, z) + + def round_to_grid(self, p): + p = self.discretise_point(p) + p_r = (p[0] * self.dx, + p[1] * self.dy, + p[2] * self.dz) + return p_r + + + def within_pml(self, p): + if (p[0] < self.pmlthickness['x0'] or + p[0] > self.nx - self.pmlthickness['xmax'] or + p[1] < self.pmlthickness['y0'] or + p[1] > self.ny - self.pmlthickness['ymax'] or + p[2] < self.pmlthickness['z0'] or + p[2] > self.nz - self.pmlthickness['zmax']): + return True + else: + return False class FDTDGrid(Grid): @@ -119,6 +143,7 @@ class FDTDGrid(Grid): self.srcsteps = [0, 0, 0] self.rxsteps = [0, 0, 0] self.snapshots = [] + self.subgrids = [] def initialise_geometry_arrays(self): """ @@ -143,6 +168,13 @@ class FDTDGrid(Grid): self.Hy = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.dtypes['float_or_double']) self.Hz = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.dtypes['float_or_double']) + + def initialise_grids(self): + """Function to call the initialisation of all grids.""" + for g in [self] + self.subgrids: + g.initialise_geometry_arrays() + g.initialise_field_arrays() + def initialise_std_update_coeff_arrays(self): """Initialise arrays for storing update coefficients.""" self.updatecoeffsE = np.zeros((len(self.materials), 5), dtype=config.dtypes['float_or_double']) @@ -362,9 +394,9 @@ def dispersion_analysis(G): minwavelength = minvelocity / results['maxfreq'] # Maximum spatial step - if '3D' in config.mode: + if '3D' in config.general['mode']: delta = max(G.dx, G.dy, G.dz) - elif '2D' in config.mode: + elif '2D' in config.general['mode']: if G.nx == 1: delta = max(G.dy, G.dz) elif G.ny == 1: diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py index bee39931..338815d7 100644 --- a/gprMax/model_build_run.py +++ b/gprMax/model_build_run.py @@ -110,6 +110,8 @@ class ModelBuildRun: else: self.reuse_geometry() + G = self.G + # Adjust position of simple sources and receivers if required if G.srcsteps[0] != 0 or G.srcsteps[1] != 0 or G.srcsteps[2] != 0: for source in itertools.chain(G.hertziandipoles, G.magneticdipoles): @@ -129,7 +131,7 @@ class ModelBuildRun: receiver.zcoord = receiver.zcoordorigin + (currentmodelrun - 1) * G.rxsteps[2] # Write files for any geometry views and geometry object outputs - if not (G.geometryviews or G.geometryobjectswrite) and args.geometry_only and config.general['messages']: + if not (G.geometryviews or G.geometryobjectswrite) and self.sim_config.geometry_only and config.general['messages']: print(Fore.RED + '\nWARNING: No geometry views or geometry objects to output found.' + Style.RESET_ALL) if config.general['messages']: print() for i, geometryview in enumerate(G.geometryviews): @@ -143,7 +145,7 @@ class ModelBuildRun: pbar.close() # If only writing geometry information - if args.geometry_only: + if self.sim_config.geometry_only: tsolve = 0 def build_geometry(self): diff --git a/gprMax/pml.py b/gprMax/pml.py index 51932dd8..d0432317 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.py @@ -445,7 +445,7 @@ def pml_information(G): pmlinfo = pmlinfo[:-2] + ' cells' return 'PML: formulation: {}, order: {}, thickness: {}'.format(G.pmlformulation, len(G.cfs), pmlinfo) -def build_pml(G, key, thickness): +def build_pml(G, key, value): """This function builds instances of the PML and calculates the initial parameters and coefficients including setting profile (based on underlying material er and mr from solid array). diff --git a/gprMax/scene.py b/gprMax/scene.py index 9b2b6596..4d9a61f6 100644 --- a/gprMax/scene.py +++ b/gprMax/scene.py @@ -5,7 +5,7 @@ from .cmds_single_use import DomainSingle from .cmds_single_use import Discretisation from .cmds_single_use import TimeWindow from .cmds_multiple import UserObjectMulti -from .subgrids.user_objects import SubGridBase +from .subgrids.user_objects import SubGridBase as SubGridUserBase from .cmds_geometry.cmds_geometry import UserObjectGeometry from .exceptions import CmdInputError from .cmds_geometry.fractal_box_builder import FractalBoxBuilder diff --git a/gprMax/user_inputs.py b/gprMax/user_inputs.py index 62ad67c0..c7a23091 100644 --- a/gprMax/user_inputs.py +++ b/gprMax/user_inputs.py @@ -66,6 +66,7 @@ class MainGridUserInput(UserInput): super().__init__(grid) def check_point(self, p, cmd_str, name=''): + """Discretise point and check its within the domain""" p = self.discretise_point(p) self.point_within_bounds(p, cmd_str, name) return p