From 7bc37d0ef5fa954509b86c461d8bb11b353b38f9 Mon Sep 17 00:00:00 2001 From: John Hartley Date: Wed, 21 Aug 2019 21:14:57 +0100 Subject: [PATCH] subgridding files --- gprMax/cmds_geometry/add_surface_roughness.py | 7 +- gprMax/{ => cmds_geometry}/build_templates.py | 0 gprMax/cmds_geometry/fractal_box.py | 5 +- gprMax/cmds_multiple.py | 4 +- gprMax/config.py | 93 ++++-- gprMax/contexts.py | 32 ++- gprMax/fields_outputs.py | 2 +- gprMax/fractals.py | 1 - gprMax/geometry_outputs.py | 2 +- gprMax/gprMax.py | 5 +- gprMax/grid.py | 25 +- gprMax/input_cmds_geometry.py | 1 + gprMax/input_cmds_multiuse.py | 18 +- gprMax/model_build_run.py | 272 ++++++++++-------- gprMax/scene.py | 5 +- gprMax/solvers.py | 36 ++- gprMax/sources.py | 1 - gprMax/subgrids/base.py | 7 +- gprMax/subgrids/factory.py | 46 --- gprMax/subgrids/hsg.py | 13 +- gprMax/subgrids/precursor_nodes.py | 11 +- gprMax/subgrids/updates.py | 200 +++++++++++++ gprMax/subgrids/user_objects.py | 53 ++-- gprMax/updates.py | 163 ++++++----- gprMax/user_inputs.py | 28 +- gprMax/utilities.py | 9 + 26 files changed, 656 insertions(+), 383 deletions(-) rename gprMax/{ => cmds_geometry}/build_templates.py (100%) delete mode 100644 gprMax/subgrids/factory.py create mode 100644 gprMax/subgrids/updates.py diff --git a/gprMax/cmds_geometry/add_surface_roughness.py b/gprMax/cmds_geometry/add_surface_roughness.py index 65289fea..227110b8 100644 --- a/gprMax/cmds_geometry/add_surface_roughness.py +++ b/gprMax/cmds_geometry/add_surface_roughness.py @@ -3,6 +3,7 @@ from .cmds_geometry import UserObjectGeometry from ..exceptions import CmdInputError from ..fractals import FractalSurface from ..utilities import round_value +import gprMax.config as config from tqdm import tqdm import numpy as np @@ -23,8 +24,8 @@ class AddSurfaceRoughness(UserObjectGeometry): p1 = self.kwargs['p1'] p2 = self.kwargs['p2'] frac_dim = self.kwargs['frac_dim'] - weighting = np.array(self.kwargs['weighting']) - limits = self.kwargs['limits'] + weighting = np.array(self.kwargs['weighting'], dtype=np.float64) + limits = np.array(self.kwargs['limits']) fractal_box_id = self.kwargs['fractal_box_id'] except KeyError: raise CmdInputError(self.__str__() + ' Incorrect parameters') @@ -122,5 +123,5 @@ class AddSurfaceRoughness(UserObjectGeometry): surface.generate_fractal_surface(grid) volume.fractalsurfaces.append(surface) - if grid.messages: + if config.general['messages']: tqdm.write('Fractal surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with fractal dimension {:g}, fractal weightings {:g}, {:g}, fractal seeding {}, and range {:g}m to {:g}m, added to {}.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, surface.dimension, surface.weighting[0], surface.weighting[1], surface.seed, limits[0], limits[1], surface.operatingonID)) diff --git a/gprMax/build_templates.py b/gprMax/cmds_geometry/build_templates.py similarity index 100% rename from gprMax/build_templates.py rename to gprMax/cmds_geometry/build_templates.py diff --git a/gprMax/cmds_geometry/fractal_box.py b/gprMax/cmds_geometry/fractal_box.py index 1d96bd07..959a37dd 100644 --- a/gprMax/cmds_geometry/fractal_box.py +++ b/gprMax/cmds_geometry/fractal_box.py @@ -2,6 +2,7 @@ from .cmds_geometry import UserObjectGeometry from ..exceptions import CmdInputError from ..fractals import FractalVolume +import gprMax.config as config from tqdm import tqdm import numpy as np @@ -43,7 +44,7 @@ class FractalBox(UserObjectGeometry): averagefractalbox = self.kwargs['averaging'] except KeyError: # if they havent specfied - go with the grid default - averagefractalbox = grid.averagevolumeobjects + averagefractalbox = False p1, p2 = uip.check_box_points(p1, p2, self.__str__()) xs, ys, zs = p1 @@ -82,7 +83,7 @@ class FractalBox(UserObjectGeometry): volume.averaging = averagefractalbox volume.mixingmodel = mixingmodel - if grid.messages: + if config.general['messages']: if volume.averaging: dielectricsmoothing = 'on' else: diff --git a/gprMax/cmds_multiple.py b/gprMax/cmds_multiple.py index aec363ea..28a16dd1 100644 --- a/gprMax/cmds_multiple.py +++ b/gprMax/cmds_multiple.py @@ -30,7 +30,7 @@ from .sources import TransmissionLine as TransmissionLineUser from .snapshots import Snapshot as SnapshotUser from .receivers import Rx as RxUser from .materials import Material as MaterialUser -from .materials import PeplinskiSoil +from .materials import PeplinskiSoil as PeplinskiSoilUser from .geometry_outputs import GeometryObjects as GeometryObjectsUser from .pml import CFSParameter from .pml import CFS @@ -838,7 +838,7 @@ class SoilPeplinski(UserObjectMulti): raise CmdInputError(self.__str__() + ' with ID {} already exists'.format(ID)) # Create a new instance of the Material class material (start index after pec & free_space) - s = PeplinskiSoil(ID, sand_fraction, clay_fraction, bulk_density, sand_density, (water_fraction_lower, water_fraction_upper)) + s = PeplinskiSoilUser(ID, sand_fraction, clay_fraction, bulk_density, sand_density, (water_fraction_lower, water_fraction_upper)) if config.general['messages']: print('Mixing model (Peplinski) used to create {} with sand fraction {:g}, clay fraction {:g}, bulk density {:g}g/cm3, sand particle density {:g}g/cm3, and water volumetric fraction {:g} to {:g} created.'.format(s.ID, s.S, s.C, s.rb, s.rs, s.mu[0], s.mu[1])) diff --git a/gprMax/config.py b/gprMax/config.py index b39075e4..618a6da0 100644 --- a/gprMax/config.py +++ b/gprMax/config.py @@ -45,8 +45,11 @@ z0 = np.sqrt(m0 / e0) # progressbars: whether to show progress bars on stdoout or not # mode: 2D TMx, 2D TMy, 2D TMz, or 3D # cpu, cuda, opencl: solver type -general = {'inputfilepath': '', 'outputfilepath': '', 'messages': True, - 'progressbars': True, 'mode': '3D', 'cpu': True, 'cuda': False, 'opencl': False} +# autotranslate: auto translate objects with main grid coordinates +# to their equivalent local grid coordinate within the subgrid. If this option is off +# users must specify sub-grid object point within the global subgrid space. +general = {'inputfilepath': 'gprMax', 'outputfilepath': 'gprMax', 'messages': True, + 'progressbars': True, 'mode': '3D', 'cpu': True, 'cuda': False, 'opencl': False, 'autotranslate': False} # Store information about host machine hostinfo = get_host_info() @@ -76,7 +79,7 @@ materials = {'maxpoles': 0, 'dispersivedtype': None, 'dispersiveCdtype': None} # which are represented as two floats # Main field arrays use floats (float_or_double) and complex numbers (complex) # Precision of float_or_double and complex: single or double for numpy and C (CUDA) arrays -precision = 'single' +precision = 'double' if precision == 'single': dtypes = {'float_or_double': np.float32, 'complex': np.complex64, @@ -92,13 +95,26 @@ class ModelConfig(): def __init__(self, sim_config, i): self.sim_config = sim_config + self.reuse_geometry = False # current model number (indexed from 0) self.i = i + parts = self.sim_config.output_file_path.parts + if not sim_config.single_model: # 1 indexed self.appendmodelnumber = str(self.i + 1) + else: + self.appendmodelnumber = '' + + # outputfilepath for specific model + self.output_file_path = Path(*parts[:-2], parts[-1] + self.appendmodelnumber) + self.output_file_path_ext = self.output_file_path.with_suffix('.out') + + # make a snapshot directory + stem = parts[-1] + '_snaps' + self.appendmodelnumber + self.snapshot_dir = Path(*parts[:-2], stem) inputfilestr_f = '\n--- Model {}/{}, input file: {}' self.inputfilestr = inputfilestr_f.format(self.i + 1, self.sim_config.model_end, self.sim_config.input_file_path) @@ -123,6 +139,7 @@ class ModelConfig(): 'current_model_run': self.i + 1, 'inputfile': self.sim_config.input_file_path.resolve()} + class SimulationConfig: def __init__(self, args): @@ -150,18 +167,28 @@ class SimulationConfig: self.mpi = args.mpi self.mpi_no_spawn = args.mpi_no_spawn self.general = {} - self.general['messages'] = True + self.general['messages'] = general['messages'] self.geometry_fixed = args.geometry_fixed self.geometry_only = args.geometry_only self.write_processed = args.write_processed + + # subgrid parameter may not exist if user uses CLI api + try: + self.subgrid = args.subgrid + except AttributeError: + # this must be CLI user. No subgrids are available + self.subgrid = False + + # scenes parameter may not exist if user uses CLI api try: self.scenes = args.scenes except AttributeError: self.scenes = [] + # set more complex parameters self.set_input_file_path() - self.set_model_start() - self.set_model_end() + self.set_output_file_path() + self.set_model_start_end() self.set_single_model() def set_single_model(self): @@ -171,26 +198,22 @@ class SimulationConfig: self.single_model = False # for example - def set_model_start(self): + def set_model_start_end(self): - # serial simulation - if not self.mpi and not self.mpi_no_spawn: - # Set range for number of models to run - if self.args.task: - # Job array feeds args.n number of single tasks - self.model_start = self.args.task - 1 - elif self.args.restart: - self.model_start = self.args.restart - 1 - else: - self.model_start = 0 + # set range for number of models to run (internally 0 index) + if self.args.task: + # Job array feeds args.n number of single tasks + modelstart = self.args.task - 1 + modelend = self.args.task + elif self.args.restart: + modelstart = self.args.restart - 1 + modelend = modelstart + self.args.n - 1 + else: + modelstart = 0 + modelend = modelstart + self.args.n - # mpi simulation - elif self.mpi: - # Set range for number of models to run - self.modelstart = self.args.restart - 1 if self.args.restart else 0 # etc... - - def set_model_end(self): - self.model_end = 1 # for now! + self.model_start = modelstart + self.model_end = modelend def set_precision(self): pass @@ -207,16 +230,32 @@ class SimulationConfig: def set_output_file_path(self): # output file can be provided by the user. if they havent provided None # use the inputfilefile path instead - if self.args.outputfile: + try: self.output_file_path = Path(self.args.outputfile) - else: + except AttributeError: self.output_file_path = Path(self.args.inputfile) + +class SimulationConfigMPI(SimulationConfig): + + def __init__(self, args): + super().__init__(args) + + def set_model_start_end(self): + # Set range for number of models to run + self.model_start = self.args.restart if self.args.restart else 1 + self.model_end = self.modelstart + self.args.n + + def create_simulation_config(args): - sc = SimulationConfig(args) + if not args.mpi and not args.mpi_no_spawn: + sc = SimulationConfig(args) + elif args.mpi: + sc = SimulationConfigMPI(args) return sc + def create_model_config(sim_config, i): mc = ModelConfig(sim_config, i) return mc diff --git a/gprMax/contexts.py b/gprMax/contexts.py index a71972e9..ca480800 100644 --- a/gprMax/contexts.py +++ b/gprMax/contexts.py @@ -20,10 +20,13 @@ from .utilities import timer from .model_build_run import ModelBuildRun import datetime from .config import create_model_config +from .solvers import create_solver +from .solvers import create_G + class Context(): - def __init__(self, sim_config, solver): + def __init__(self, sim_config): """Context for the model to run in. Sub-class this with contexts i.e. an MPI context. @@ -32,7 +35,6 @@ class Context(): solver (Solver): FDTD general solver object. """ self.sim_config = sim_config - self.solver = solver self.model_range = range(sim_config.model_start, sim_config.model_end) self.tsimend = 0 @@ -58,13 +60,26 @@ class NoMPIContext(Context): def _run(self): """Specialise how the models are farmed out.""" + for i in self.model_range: - # create the model configuration model_config = create_model_config(self.sim_config, i) - model = ModelBuildRun(self.solver, self.sim_config, model_config) + # always create a solver for the first model + # the next model to run only gets a new solver if the + # geometry is not re used. + if i != 0 and self.sim_config.geometry_fixed: + model_config.reuse_geometry = True + else: + G = create_G(self.sim_config) + + model = ModelBuildRun(G, self.sim_config, model_config) model.build() + solver = create_solver(G, self.sim_config) + + if not self.sim_config.geometry_only: + model.run_model(solver) + def make_time_report(self): """Function to specialise the time reporting for the standard Simulation context.""" @@ -82,6 +97,7 @@ class MPIContext(Context): def make_time_report(self): pass + class MPINoSpawnContext(Context): def _run(self): @@ -91,13 +107,13 @@ class MPINoSpawnContext(Context): pass -def create_context(sim_config, solver): +def create_context(sim_config): """Create a context in which to run the simulation. i.e MPI.""" if sim_config.mpi_no_spawn: - context = MPIContext(sim_config, solver) + context = MPIContext(sim_config) elif sim_config.mpi: - context = MPINoSpawnContext(sim_config, solver) + context = MPINoSpawnContext(sim_config) else: - context = NoMPIContext(sim_config, solver) + context = NoMPIContext(sim_config) return context diff --git a/gprMax/fields_outputs.py b/gprMax/fields_outputs.py index 57b8f0b1..7a9bdd39 100644 --- a/gprMax/fields_outputs.py +++ b/gprMax/fields_outputs.py @@ -92,7 +92,7 @@ __global__ void store_outputs(int NRX, int iteration, const int* __restrict__ rx """) -def write_hdf5_outputfile(outputfile, Ex, Ey, Ez, Hx, Hy, Hz, G): +def write_hdf5_outputfile(outputfile, G): """Write an output file in HDF5 format. Args: diff --git a/gprMax/fractals.py b/gprMax/fractals.py index 3822e016..f2ff73c2 100644 --- a/gprMax/fractals.py +++ b/gprMax/fractals.py @@ -95,7 +95,6 @@ class FractalSurface(object): # Generate fractal generate_fractal2D(surfacedims[0], surfacedims[1], config.hostinfo['ompthreads'], self.b, self.weighting, v1, A, self.fractalsurface) - # Shift the zero frequency component to start of the array self.fractalsurface = fftpack.ifftshift(self.fractalsurface) # Take the real part (numerical errors can give rise to an imaginary part) of the IFFT diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py index fe8da4ea..ee0139bb 100644 --- a/gprMax/geometry_outputs.py +++ b/gprMax/geometry_outputs.py @@ -39,7 +39,7 @@ class GeometryView(object): else: byteorder = 'BigEndian' - def __init__(self, xs=None, ys=None, zs=None, xf=None, yf=None, zf=None, dx=None, dy=None, dz=None, filename=None, fileext=None): + def __init__(self, xs=None, ys=None, zs=None, xf=None, yf=None, zf=None, dx=None, dy=None, dz=None, filename=None, fileext=None, grid=None): """ Args: xs, xf, ys, yf, zs, zf (int): Extent of the volume in cells. diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index 09405e11..3830171b 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.py @@ -34,6 +34,7 @@ def api( mpi_no_spawn=False, mpicomm=None, gpu=None, + subgrid=None, benchmark=False, geometry_only=False, geometry_fixed=False, @@ -56,6 +57,7 @@ def api( args.mpi_no_spawn = mpi_no_spawn args.mpicomm = mpicomm args.gpu = gpu + args.subgrid=subgrid args.benchmark = benchmark args.geometry_only = geometry_only args.geometry_fixed = geometry_fixed @@ -88,6 +90,5 @@ def main(): def run_main(args): sim_config = create_simulation_config(args) - solver = create_solver(sim_config) - context = create_context(sim_config, solver) + context = create_context(sim_config) context.run() diff --git a/gprMax/grid.py b/gprMax/grid.py index edfc1b55..6c9ceca8 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -33,6 +33,10 @@ from .utilities import fft_power from .utilities import human_size from .utilities import round_value +import decimal as d +from scipy.constants import c + + class Grid(object): """Generic grid/mesh.""" @@ -145,6 +149,8 @@ class FDTDGrid(Grid): self.snapshots = [] self.subgrids = [] self.gpu = None + self.name = 'Main' + self.outputdirectory = '' def initialise_geometry_arrays(self): """ @@ -296,12 +302,23 @@ class FDTDGrid(Grid): def reset_fields(self): # Clear arrays for field components - G.initialise_field_arrays() + self.initialise_field_arrays() # Clear arrays for fields in PML - for pml in G.pmls: + for pml in self.pmls: pml.initialise_field_arrays() + def calculate_dt(self): + self.dt = 1 / (c * np.sqrt((1 / self.dx) * (1 / self.dx) + (1 / self.dy) * (1 / self.dy) + (1 / self.dz) * (1 / self.dz))) + + def round_time_step(self): + # 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. + # 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. + self.dt = round_value(self.dt, decimalplaces=d.getcontext().prec - 1) + + + def dispersion_analysis(G): """ Analysis of numerical dispersion (Taflove et al, 2005, p112) - @@ -477,3 +494,7 @@ def Iz(x, y, z, Hx, Hy, Hz, G): Iz = G.dx * (Hx[x, y - 1, z] - Hx[x, y, z]) + G.dy * (Hy[x, y, z] - Hy[x - 1, y, z]) return Iz + + +class GPUGrid(Grid): + pass diff --git a/gprMax/input_cmds_geometry.py b/gprMax/input_cmds_geometry.py index 17ffb18d..8010f223 100644 --- a/gprMax/input_cmds_geometry.py +++ b/gprMax/input_cmds_geometry.py @@ -36,6 +36,7 @@ from .utilities import round_value from tqdm import tqdm +import numpy as np def process_geometrycmds(geometry): diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py index 27072274..23aaac9d 100644 --- a/gprMax/input_cmds_multiuse.py +++ b/gprMax/input_cmds_multiuse.py @@ -43,7 +43,7 @@ from .cmds_multiple import Snapshot from .cmds_multiple import AddDebyeDispersion from .cmds_multiple import AddLorentzDispersion from .cmds_multiple import AddDrudeDispersion -from .cmds_multiple import PeplinskiSoil +from .cmds_multiple import SoilPeplinski from .cmds_multiple import GeometryView from .cmds_multiple import GeometryObjectsWrite from .cmds_multiple import PMLCFS @@ -268,15 +268,13 @@ def process_multicmds(multicmds): if len(tmp) != 7: raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires at exactly seven parameters') - - soil = PeplinskiSoil( - sand_fraction=float(tmp[0]), - clay_fraction=float(tmp[1]), - bulk_density=float(tmp[2]), - sand_density=float(tmp[3]), - water_fraction_lower=float(tmp[4]), - water_fraction_upper=float(tmp[5]), - ID=tmp[6]) + soil = SoilPeplinski(sand_fraction=float(tmp[0]), + clay_fraction=float(tmp[1]), + bulk_density=float(tmp[2]), + sand_density=float(tmp[3]), + water_fraction_lower=float(tmp[4]), + water_fraction_upper=float(tmp[5]), + ID=tmp[6]) scene_objects.append(soil) # Geometry views (creates VTK-based geometry files) diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py index 4cab2cf1..41f205e9 100644 --- a/gprMax/model_build_run.py +++ b/gprMax/model_build_run.py @@ -31,6 +31,7 @@ import cython import numpy as np from terminaltables import SingleTable from tqdm import tqdm +from pathlib import Path import gprMax.config as config from .cuda.fields_updates import kernel_template_fields @@ -39,7 +40,6 @@ from .cuda.source_updates import kernel_template_sources from .cython.yee_cell_build import build_electric_components from .cython.yee_cell_build import build_magnetic_components from .exceptions import GeneralError -from .fields_outputs import store_outputs from .fields_outputs import kernel_template_store_outputs from .fields_outputs import write_hdf5_outputfile from .grid import FDTDGrid @@ -69,24 +69,21 @@ from .utilities import human_size from .utilities import open_path_file from .utilities import round32 from .utilities import timer +from .utilities import Printer from .scene import Scene +from .solvers import create_solver -class Printer(): - def __init__(self, sim_config): - self.printing = sim_config.general['messages'] - - def print(self, str): - if self.printing: - print(str) class ModelBuildRun: - def __init__(self, solver, sim_config, model_config): - self.solver = solver + def __init__(self, G, sim_config, model_config): + self.G = G self.sim_config = sim_config self.model_config = model_config - self.G = solver.get_G() + self.printer = Printer(sim_config) + # Monitor memory usage + self.p = None def build(self): """Runs a model - processes the input file; builds the Yee cells; calculates update coefficients; runs main FDTD loop. @@ -104,13 +101,13 @@ class ModelBuildRun: tsolve (int): Length of time (seconds) of main FDTD calculations """ # Monitor memory usage - p = psutil.Process() + self.p = psutil.Process() # Normal model reading/building process; bypassed if geometry information to be reused - if not self.sim_config.geometry_fixed: - self.build_geometry() - else: + if self.model_config.reuse_geometry: self.reuse_geometry() + else: + self.build_geometry() G = self.G @@ -137,7 +134,7 @@ class ModelBuildRun: 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): - geometryview.set_filename(appendmodelnumber) + geometryview.set_filename(self.model_config.appendmodelnumber) pbar = tqdm(total=geometryview.datawritesize, unit='byte', unit_scale=True, desc='Writing geometry view file {}/{}, {}'.format(i + 1, len(G.geometryviews), os.path.split(geometryview.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) geometryview.write_vtk(G, pbar) pbar.close() @@ -153,6 +150,7 @@ class ModelBuildRun: def build_geometry(self): model_config = self.model_config sim_config = self.sim_config + G = self.G printer = Printer(sim_config) @@ -160,19 +158,17 @@ class ModelBuildRun: scene = self.build_scene() - # print PML information - printer.print(pml_information(G)) + # combine available grids + grids = [G] + G.subgrids + gridbuilders = [GridBuilder(grid, self.printer) for grid in grids] - self.build_pmls() - self.build_components() - - # update grid for tm modes - self.tm_grid_update() - - self.update_voltage_source_materials() - - # Initialise arrays of update coefficients to pass to update functions - G.initialise_std_update_coeff_arrays() + for gb in gridbuilders: + gb.printer.print(pml_information(gb.grid)) + gb.build_pmls() + gb.build_components() + gb.tm_grid_update() + gb.update_voltage_source_materials() + gb.grid.initialise_std_update_coeff_arrays() # Set datatype for dispersive arrays if there are any dispersive materials. if config.materials['maxpoles'] != 0: @@ -189,7 +185,8 @@ class ModelBuildRun: G.memory_check() printer.print('\nMemory (RAM) required - updated (dispersive): ~{}\n'.format(human_size(G.memoryusage))) - G.initialise_dispersive_arrays(config.materials['dispersivedtype']) + for gb in gridbuilders: + gb.grid.initialise_dispersive_arrays(config.materials['dispersivedtype']) # Check there is sufficient memory to store any snapshots if G.snapshots: @@ -202,15 +199,8 @@ class ModelBuildRun: printer.print('\nMemory (RAM) required - updated (snapshots): ~{}\n'.format(human_size(G.memoryusage))) - # Process complete list of materials - calculate update coefficients, - # store in arrays, and build text list of materials/properties - materialsdata = process_materials(G) - materialstable = SingleTable(materialsdata) - materialstable.outer_border = False - materialstable.justify_columns[0] = 'right' - - printer.print('\nMaterials:') - printer.print(materialstable.table) + for gb in gridbuilders: + gb.build_materials() # Check to see if numerical dispersion might be a problem results = dispersion_analysis(G) @@ -223,25 +213,16 @@ class ModelBuildRun: 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'])) - # set the dispersive update functions based on the model configuration - props = self.solver.updates.adapt_dispersive_config(config) - self.solver.updates.set_dispersive_updates(props) def reuse_geometry(self): - printer = Printer(model_config) - inputfilestr = '\n--- Model {}/{}, input file (not re-processed, i.e. geometry fixed): {}'.format(currentmodelrun, modelend, model_config.input_file_path) - printer.print(Fore.GREEN + '{} {}\n'.format(model_config.inputfilestr, '-' * (get_terminal_width() - 1 - len(model_config.inputfilestr))) + Style.RESET_ALL) - self.G.reset_fields() - - def tm_grid_update(self): - if '2D TMx' == config.general['mode']: - self.G.tmx() - elif '2D TMy' == config.general['mode']: - self.G.tmy() - elif '2D TMz' == config.general['mode']: - self.G.tmz() + G = self.G + inputfilestr = '\n--- Model {}/{}, input file (not re-processed, i.e. geometry fixed): {}'.format(self.model_config.appendmodelnumber, self.sim_config.model_end, self.sim_config.input_file_path) + self.printer.print(Fore.GREEN + '{} {}\n'.format(self.model_config.inputfilestr, '-' * (get_terminal_width() - 1 - len(inputfilestr))) + Style.RESET_ALL) + for grid in [G] + G.subgrids: + grid.reset_fields() def build_scene(self): + G = self.G # api for multiple scenes / model runs scene = self.model_config.get_scene() @@ -249,96 +230,135 @@ class ModelBuildRun: if not scene: scene = Scene() # parse the input file into user objects and add them to the scene - scene = parse_hash_commands(self.model_config, self.G, scene) + scene = parse_hash_commands(self.model_config, G, scene) # Creates the internal simulation objects. - scene.create_internal_objects(self.G) + scene.create_internal_objects(G) return scene - def build_pmls(self): - # build the PMLS - pbar = tqdm(total=sum(1 for value in self.G.pmlthickness.values() if value > 0), desc='Building PML boundaries', ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) + def create_output_directory(self): - for pml_id, thickness in self.G.pmlthickness.items(): - build_pml(self.G, pml_id, thickness) + if self.G.outputdirectory: + # Check and set output directory and filename + try: + os.mkdir(self.G.outputdirectory) + self.printer.print('\nCreated output directory: {}'.format(self.G.outputdirectory)) + except FileExistsError: + pass + # modify the output path (hack) + self.model_config.output_file_path_ext = Path(self.G.outputdirectory, self.model_config.output_file_path_ext) + + + def run_model(self, solver): + + G = self.G + + self.create_output_directory() + self.printer.print('\nOutput file: {}\n'.format(self.model_config.output_file_path_ext)) + + tsolve = self.solve(solver) + + # Write an output file in HDF5 format + write_hdf5_outputfile(self.model_config.output_file_path_ext, G) + + # Write any snapshots to file + if G.snapshots: + # Create directory and construct filename from user-supplied name and model run number + snapshotdir = self.model_config.snapshot_dir + if not os.path.exists(snapshotdir): + os.mkdir(snapshotdir) + + self.printer.print() + for i, snap in enumerate(G.snapshots): + snap.filename = snapshotdir + snap.basefilename + '.vti' + pbar = tqdm(total=snap.vtkdatawritesize, leave=True, unit='byte', unit_scale=True, desc='Writing snapshot file {} of {}, {}'.format(i + 1, len(G.snapshots), os.path.split(snap.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) + snap.write_vtk_imagedata(pbar, G) + pbar.close() + self.printer.print() + + memGPU = '' + if config.cuda['gpus']: + memGPU = ' host + ~{} GPU'.format(human_size(self.solver.get_memsolve())) + + self.printer.print('\nMemory (RAM) used: ~{}{}'.format(human_size(self.p.memory_full_info().uss), memGPU)) + self.printer.print('Solving time [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=tsolve))) + + return tsolve + + def solve(self, solver): + """ + Solving using FDTD method on CPU. Parallelised using Cython (OpenMP) for + electric and magnetic field updates, and PML updates. + + Args: + currentmodelrun (int): Current model run number. + modelend (int): Number of last model to run. + G (class): Grid class instance - holds essential parameters describing the model. + + Returns: + tsolve (float): Time taken to execute solving (seconds) + """ + G = self.G + + if config.general['messages']: + iterator = tqdm(range(G.iterations), desc='Running simulation, model ' + str(self.model_config + .i + 1) + '/' + str(self.sim_config.model_end), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) + else: + iterator = range(0, G.iterations) + + tsolve = solver.solve(iterator) + + return tsolve + + +class GridBuilder: + def __init__(self, grid, printer): + self.grid = grid + self.printer = printer + + def build_pmls(self): + + grid = self.grid + # build the PMLS + pbar = tqdm(total=sum(1 for value in grid.pmlthickness.values() if value > 0), desc='Building {} Grid PML boundaries'.format(self.grid.name), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) + + for pml_id, thickness in grid.pmlthickness.items(): + build_pml(grid, pml_id, thickness) pbar.update() pbar.close() def build_components(self): # Build the model, i.e. set the material properties (ID) for every edge # of every Yee cell - G = self.G - printer = Printer(self.sim_config) - printer.print('') - pbar = tqdm(total=2, desc='Building main grid', ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) - build_electric_components(G.solid, G.rigidE, G.ID, G) + self.printer.print('') + pbar = tqdm(total=2, desc='Building {} grid'.format(self.grid.name), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) + build_electric_components(self.grid.solid, self.grid.rigidE, self.grid.ID, self.grid) pbar.update() - build_magnetic_components(G.solid, G.rigidH, G.ID, G) + build_magnetic_components(self.grid.solid, self.grid.rigidH, self.grid.ID, self.grid) pbar.update() pbar.close() + def tm_grid_update(self): + if '2D TMx' == config.general['mode']: + self.grid.tmx() + elif '2D TMy' == config.general['mode']: + self.grid.tmy() + elif '2D TMz' == config.general['mode']: + self.grid.tmz() + def update_voltage_source_materials(self): # Process any voltage sources (that have resistance) to create a new # material at the source location - for voltagesource in self.G.voltagesources: - voltagesource.create_material(self.G) + for voltagesource in self.grid.voltagesources: + voltagesource.create_material(self.grid) + def build_materials(self): + # Process complete list of materials - calculate update coefficients, + # store in arrays, and build text list of materials/properties + materialsdata = process_materials(self.grid) + materialstable = SingleTable(materialsdata) + materialstable.outer_border = False + materialstable.justify_columns[0] = 'right' - def run_model(self): - - # Check and set output directory and filename - if not os.path.isdir(config.outputfilepath): - os.mkdir(config.outputfilepath) - printer.print('\nCreated output directory: {}'.format(config.outputfilepath)) - - outputfile = os.path.join(config.outputfilepath, os.path.splitext(os.path.split(config.inputfilepath)[1])[0] + appendmodelnumber + '.out') - printer.print('\nOutput file: {}\n'.format(outputfile)) - - tsolve, memsolve = solve(currentmodelrun, modelend, G) - - # Write an output file in HDF5 format - write_hdf5_outputfile(outputfile, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G) - - # Write any snapshots to file - if G.snapshots: - # Create directory and construct filename from user-supplied name and model run number - snapshotdir = os.path.splitext(config.inputfilepath)[0] + '_snaps' + appendmodelnumber - if not os.path.exists(snapshotdir): - os.mkdir(snapshotdir) - - if config.general['messages']: print() - for i, snap in enumerate(G.snapshots): - snap.filename = os.path.abspath(os.path.join(snapshotdir, snap.basefilename + '.vti')) - pbar = tqdm(total=snap.vtkdatawritesize, leave=True, unit='byte', unit_scale=True, desc='Writing snapshot file {} of {}, {}'.format(i + 1, len(G.snapshots), os.path.split(snap.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) - snap.write_vtk_imagedata(pbar, G) - pbar.close() - if config.general['messages']: print() - - memGPU = '' - if config.cuda['gpus']: - memGPU = ' host + ~{} GPU'.format(human_size(memsolve)) - - printer.print('\nMemory (RAM) used: ~{}{}'.format(human_size(p.memory_full_info().uss), memGPU)) - printer.print('Solving time [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=tsolve))) - - - return tsolve - -def solve(solver, sim_config, model_config): - """ - Solving using FDTD method on CPU. Parallelised using Cython (OpenMP) for - electric and magnetic field updates, and PML updates. - - Args: - currentmodelrun (int): Current model run number. - modelend (int): Number of last model to run. - G (class): Grid class instance - holds essential parameters describing the model. - - Returns: - tsolve (float): Time taken to execute solving (seconds) - """ - tsolvestart = timer() - - solver.solve() - - return tsolve + self.printer.print('\n{} Grid Materials:'.format(self.grid.name)) + self.printer.print(materialstable.table) diff --git a/gprMax/scene.py b/gprMax/scene.py index 4d9a61f6..a298b887 100644 --- a/gprMax/scene.py +++ b/gprMax/scene.py @@ -47,6 +47,7 @@ class Scene: else: return False + # subgrid user objects subgrid_cmds = list(filter(func, self.multiple_cmds)) # iterate through the user command objects under the subgrid user object @@ -55,8 +56,8 @@ class Scene: # object. this reference allows the multi and geo user objects # to build in the correct subgrid. sg = sg_cmd.subgrid - self.process_cmds(sg_cmd.multiple_cmds, sg) - self.process_cmds(sg_cmd.geometry_cmds, sg) + self.process_cmds(sg_cmd.children_multiple, sg) + self.process_cmds(sg_cmd.children_geometry, sg) def process_cmds(self, commands, grid): cmds_sorted = sorted(commands, key=lambda cmd: cmd.order) diff --git a/gprMax/solvers.py b/gprMax/solvers.py index fa978c03..2f7d9024 100644 --- a/gprMax/solvers.py +++ b/gprMax/solvers.py @@ -18,21 +18,41 @@ from gprMax.updates import CPUUpdates from gprMax.updates import GPUUpdates from gprMax.utilities import timer +from .grid import FDTDGrid +from .grid import GPUGrid +import gprMax.config as config +from .subgrids.updates import create_subgrid_updates -def create_solver(sim_config): +def create_G(sim_config): """Returns the configured solver.""" if sim_config.gpu: - from .grid import GPUGrid G = GPUGrid() - updates = GPUUpdates(G) - else: - from .grid import FDTDGrid + elif sim_config.subgrid: G = FDTDGrid() + else: + G = FDTDGrid() + + return G + + +def create_solver(G, sim_config): + """Returns the configured solver.""" + if sim_config.gpu: + updates = GPUUpdates(G) + elif sim_config.subgrid: + updates = create_subgrid_updates(G) + else: updates = CPUUpdates(G) solver = Solver(updates) + # a large range of function exist to advance the time step for dispersive + # materials. The correct function is set here based on the + # the required numerical precision and dispersive material type. + props = updates.adapt_dispersive_config(config) + updates.set_dispersive_updates(props) + return solver @@ -49,17 +69,15 @@ class Solver: iterator (iterator): can be range() or tqdm() """ self.updates = updates - #self.iterator = iterator def get_G(self): return self.updates.G - def solve(self): + def solve(self, iterator): """Time step the FDTD model.""" tsolvestart = timer() - for iteration in self.iterator: - + for iteration in iterator: self.updates.store_outputs(iteration) self.updates.store_snapshots(iteration) self.updates.update_magnetic() diff --git a/gprMax/sources.py b/gprMax/sources.py index 5665ffa0..4ce3c26b 100644 --- a/gprMax/sources.py +++ b/gprMax/sources.py @@ -49,7 +49,6 @@ class Source(object): Args: G (class): Grid class instance - 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']) diff --git a/gprMax/subgrids/base.py b/gprMax/subgrids/base.py index 4bbaef84..4bb0298e 100644 --- a/gprMax/subgrids/base.py +++ b/gprMax/subgrids/base.py @@ -17,8 +17,6 @@ # along with gprMax. If not, see . from ..grid import FDTDGrid - -from scipy.constants import c import numpy as np @@ -63,9 +61,8 @@ class SubGridBase(FDTDGrid): self.n_boundary_cells_z = d_to_pml + self.pmlthickness['z0'] # vectorise coordinates - self.p0 = np.array(self.i0, self.j0, self.k0) - self.n_boundary_cells_p = np.array(self.n_boundary_cells_x, - self.n_boundary_cells_y, self.n_boundary_cells_z) + #self.p0 = np.array(self.i0, self.j0, self.k0) + self.n_boundary_cells_p = np.array([self.n_boundary_cells_x, self.n_boundary_cells_y, self.n_boundary_cells_z]) # interpolation scheme self.interpolation = kwargs['interpolation'] diff --git a/gprMax/subgrids/factory.py b/gprMax/subgrids/factory.py deleted file mode 100644 index a05963c8..00000000 --- a/gprMax/subgrids/factory.py +++ /dev/null @@ -1,46 +0,0 @@ -# Copyright (C) 2015-2019: The University of Edinburgh -# Authors: Craig Warren and Antonis Giannopoulos -# -# This file is part of gprMax. -# -# gprMax is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# gprMax is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with gprMax. If not, see . - -from .subgrid_solver import SubGridSolver -from .cpu_solvers import CPUSolver -from .subgrid_updaters import SubgridUpdater -from ..config import get_iterations -from ..config import filter -from precursor_nodes import PrecursorNodes -from precursor_nodes import PrecursorNodesFiltered - - -def create_solver(G, config): - """Factory method to return a solver instance for a CPU or subgrid simulation""" - - iterations = get_iterations(config, G) - - if filter: - from precursor_nodes import PrecursorNodesFiltered as PrecursorNodes - else: - from precursor_nodes import PrecursorNodes - - updaters = [] - for sg in G.subgrids: - precursors = PrecursorNodes(G, sg) - sgu = SubgridUpdater(sg, precursors, G) - updaters.append(sgu) - - solver = SubGridSolver(G, updaters, iterations) - - return solver diff --git a/gprMax/subgrids/hsg.py b/gprMax/subgrids/hsg.py index 87446b43..3327b40a 100644 --- a/gprMax/subgrids/hsg.py +++ b/gprMax/subgrids/hsg.py @@ -34,14 +34,6 @@ class SubGridHSG(SubGridBase): super().__init__(**kwargs) self.gridtype = SubGridHSG.gridtype - # upper and lower indices for the OS. - self.i_l = self.i0 - self.is_os_sep - self.i_u = self.i1 + self.is_os_sep - self.j_l = self.j0 - self.is_os_sep - self.j_u = self.j1 + self.is_os_sep - self.k_l = self.k0 - self.is_os_sep - self.k_u = self.k1 + self.is_os_sep - def memory_usage(self): pass @@ -127,13 +119,16 @@ class SubGridHSG(SubGridBase): cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, self.i_l, self.i_u + 1, self.j_l, self.j_u, self.k_l - 1, self.k_u, self.nwz, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ey, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads) def __str__(self): + + self.memory_estimate_basic() + s = '\n' s += Fore.CYAN s += 'Sub Grid HSG\n' s += 'Name: {}\n'.format(self.name) s += 'dx, dy, dz: {}m {}m {}m\n'.format(self.dx, self.dy, self.dz) s += 'dt: {}s\n'.format(self.dt) - s += 'Memory Estimate: {}\n'.format(human_size(self.memory_usage())) + s += 'Memory Estimate: {}\n'.format(human_size(self.memoryusage)) s += 'Position: ({}m, {}m, {}m), ({}m, {}m, {}m)\n'.format(self.x1, self.y1, self.z1, diff --git a/gprMax/subgrids/precursor_nodes.py b/gprMax/subgrids/precursor_nodes.py index ac37e5db..a57686ce 100644 --- a/gprMax/subgrids/precursor_nodes.py +++ b/gprMax/subgrids/precursor_nodes.py @@ -16,10 +16,8 @@ # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -from .config import sg_interp_d import numpy as np from scipy import interpolate -import sys def calculate_weighting_coefficients(x1, x): @@ -37,7 +35,7 @@ class PrecusorNodesBase: self.nwy = sub_grid.nwy self.nwz = sub_grid.nwz self.sub_grid = sub_grid - self.interpolation = sg_interp_d + self.interpolation = sub_grid.interpolation self.Hx = fdtd_grid.Hx self.Hy = fdtd_grid.Hy @@ -242,7 +240,7 @@ class PrecusorNodesBase: # Grab the main grid fields used to interpolate across the IS # f = self.Hi[slice] - f1, f2 = f1_f2_magnetic(obj) + f_1, f_2 = self.f1_f2_magnetic(obj) w = obj[-2] c1, c2 = calculate_weighting_coefficients(w, self.ratio) # transverse interpolated h field @@ -262,7 +260,7 @@ class PrecusorNodesBase: self.update_previous_timestep_fields(self.fn_e) for obj in self.electric_slices: - f_m = self.f_m(obj) + f_m = self.f_m_electric(obj) f_i = self.interpolate_to_sub_grid(f_m, obj[1]) # discard the outer nodes only required for interpolation f = f_i[self.ratio:-self.ratio, self.ratio:-self.ratio] @@ -277,7 +275,6 @@ class PrecursorNodesFiltered(PrecusorNodesBase): self.initialize_magnetic_slices_array() self.initialize_electric_slices_array() - def initialize_magnetic_slices_array(self): # Array contains the indices at which the main grid should be sliced @@ -620,7 +617,7 @@ class PrecursorNodes(PrecusorNodesBase): def f1_f2_magnetic(self, obj): f_1 = obj[-1][obj[2]] f_2 = obj[-1][obj[3]] - return f1, f2 + return f_1, f_2 def f_m_electric(self, obj): return obj[-1][obj[2]] diff --git a/gprMax/subgrids/updates.py b/gprMax/subgrids/updates.py new file mode 100644 index 00000000..2fb957f4 --- /dev/null +++ b/gprMax/subgrids/updates.py @@ -0,0 +1,200 @@ +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . +from ..updates import CPUUpdates + + +def create_subgrid_updates(G): + # collection responsible for updating each subgrid + updaters = [] + + # precursor nodes can use low pass filtering at the IS. Filters are either + # on or off for all subgrids + if G.subgrids[0].filter: + from .precursor_nodes import PrecursorNodesFiltered as PrecursorNodes + else: + from .precursor_nodes import PrecursorNodes + + # make an updater and precursor nodes for each sub grid + for sg in G.subgrids: + precursors = PrecursorNodes(G, sg) + sgu = SubgridUpdater(sg, precursors, G) + updaters.append(sgu) + + # responsible for updating the subgridding simulation as a whole + updates = SubgridUpdates(G, updaters) + + return updates + + +class SubgridUpdates(CPUUpdates): + """Top level subgrid updater. Manages the collection of subgrids.""" + + def __init__(self, G, subgrid_updaters): + super().__init__(G) + self.subgrid_updaters = subgrid_updaters + + def hsg_1(self): + """Method to update the subgrids over the first phase""" + for sg_updater in self.subgrid_updaters: + sg_updater.hsg_1() + + def hsg_2(self): + """Method to update the subgrids over the second phase""" + for sg_updater in self.subgrid_updaters: + sg_updater.hsg_2() + + def update_electric_a(self): + super().update_electric_a() + self.hsg_1() + + def update_magnetic(self): + super().update_magnetic() + self.hsg_2() + + def set_dispersive_updates(self, props): + # set the dispersive update functions for the main grid updates + super().set_dispersive_updates(props) + # set the same dispersive update functions for all fields in the subgrids + for updater in self.subgrid_updaters: + updater.set_dispersive_updates(props) + + +class SubgridUpdater(CPUUpdates): + """Class to handle updating the electric and magnetic fields of an HSG + subgrid. The IS, OS, subgrid region and the electric/magnetic sources are updated + using the precursor regions. + """ + + def __init__(self, subgrid, precursors, G): + """ + Args: + subgrid (SubGrid3d): Subgrid to be updated + precursors (PrecursorNodes): Precursor nodes associated with + the subgrid + G (class): Grid class instance - holds essential parameters + describing the model. + """ + super().__init__(subgrid) + self.precursors = precursors + self.source_iteration = 0 + self.G = G + + def hsg_1(self): + """This is the first half of the subgrid update. Takes the time step + up to the main grid magnetic update""" + precursors = self.precursors + + precursors.update_electric() + + # the subgrid + sub_grid = self.grid + + upper_m = int(sub_grid.ratio / 2 - 0.5) + + for m in range(1, upper_m + 1): + + # store_outputs(self.grid) + # STD update, interpolate inc. field in time, apply correction + self.update_electric_a() + + for pml in sub_grid.pmls: + pml.update_electric(sub_grid) + self.update_sub_grid_electric_sources() + precursors.interpolate_magnetic_in_time(int(m + sub_grid.ratio / 2 - 0.5)) + sub_grid.update_electric_is(precursors) + self.update_electric_b() + + # STD update, interpolate inc. field in time, apply correction + self.update_magnetic() + + for pml in sub_grid.pmls: + pml.update_magnetic(sub_grid) + precursors.interpolate_electric_in_time(m) + sub_grid.update_magnetic_is(precursors) + self.update_sub_grid_magnetic_sources() + + # store_outputs(self.grid) + self.update_electric_a() + for pml in sub_grid.pmls: + pml.update_electric(sub_grid) + + self.update_sub_grid_electric_sources() + + precursors.calc_exact_magnetic_in_time() + sub_grid.update_electric_is(precursors) + self.update_electric_b() + + sub_grid.update_electric_os(self.G) + + def hsg_2(self): + """This is the first half of the subgrid update. Takes the time step + up to the main grid electric update""" + sub_grid = self.grid + precursors = self.precursors + + precursors.update_magnetic() + + upper_m = int(sub_grid.ratio / 2 - 0.5) + + for m in range(1, upper_m + 1): + + self.update_magnetic() + + for pml in sub_grid.pmls: + pml.update_magnetic(sub_grid) + precursors.interpolate_electric_in_time(int(m + sub_grid.ratio / 2 - 0.5)) + sub_grid.update_magnetic_is(precursors) + + self.update_sub_grid_magnetic_sources() + + # store_outputs(self.grid) + self.update_electric_a() + + for pml in sub_grid.pmls: + pml.update_electric(sub_grid) + self.update_sub_grid_electric_sources() + + precursors.interpolate_magnetic_in_time(m) + sub_grid.update_electric_is(precursors) + self.update_electric_b() + + self.update_magnetic() + + for pml in sub_grid.pmls: + pml.update_magnetic(sub_grid) + + precursors.calc_exact_electric_in_time() + sub_grid.update_magnetic_is(precursors) + self.update_sub_grid_magnetic_sources() + + sub_grid.update_magnetic_os(self.G) + + def update_sub_grid_electric_sources(self): + """Update any electric sources in the subgrid""" + sg = self.grid + for source in sg.voltagesources + sg.transmissionlines + sg.hertziandipoles: + source.update_electric(self.source_iteration, sg.updatecoeffsE, sg.ID, sg.Ex, sg.Ey, sg.Ez, sg) + + self.source_iteration += 1 + + def update_sub_grid_magnetic_sources(self): + """Update any magnetic sources in the subgrid""" + sg = self.grid + for source in sg.transmissionlines + sg.magneticdipoles: + source.update_magnetic(self.source_iteration, sg.updatecoeffsH, sg.ID, + sg.Hx, sg.Hy, sg.Hz, sg) diff --git a/gprMax/subgrids/user_objects.py b/gprMax/subgrids/user_objects.py index 36793e3f..624115c2 100644 --- a/gprMax/subgrids/user_objects.py +++ b/gprMax/subgrids/user_objects.py @@ -4,6 +4,7 @@ from ..exceptions import CmdInputError from ..cmds_multiple import UserObjectMulti from ..cmds_geometry.cmds_geometry import UserObjectGeometry from ..cmds_multiple import Rx +import gprMax.config as config from copy import copy @@ -38,28 +39,23 @@ class SubGridBase(UserObjectMulti): def set_discretisation(self, sg, grid): """Set the spatial discretisation.""" - sg.dx = grid.dx / sg.ratio - sg.dy = grid.dy / sg.ratio - sg.dz = grid.dz / sg.ratio + sg.dl = grid.dl / sg.ratio + sg.dx, sg.dy, sg.dz = grid.dl def set_main_grid_indices(self, sg, grid, uip, p1, p2): """Set subgrid indices related to main grid placement.""" - # Main grid indices of the sub grid. These are dummy indices. They are - # not user internal except for printing to the user - sg.i0_u, sg.j0_u, sg.k0_u = p1 - sg.i1_u, sg.j1_u, sg.k1_u = p2 - # The actual sub gridded area (IS index) is 4 cells in - sg.i0, sg.j0, sg.k0 = np.add([sg.i0_u, sg.j0_u, sg.k0_u], sg.is_os_sep) - sg.i1, sg.j1, sg.k1 = np.subtract([sg.i1_u, sg.j1_u, sg.k1_u], sg.is_os_sep) + # IS indices + sg.i0, sg.j0, sg.k0 = p1 + sg.i1, sg.j1, sg.k1 = p2 - # Main grid indices of the sub grid. These are dummy indices. They are - # not user internal except for printing to the user - sg.x1_u, sg.y1_u, sg.z1_u = uip.round_to_grid(p1) - sg.x2_u, sg.y2_u, sg.z2_u = uip.round_to_grid(p2) + # OS indices + sg.i_l, sg.j_l, sg.k_l = p1 - sg.is_os_sep + sg.i_u, sg.j_u, sg.k_u = p2 + sg.is_os_sep - sg.x1, sg.y1, sg.z1 = np.add([sg.x1_u, sg.y1_u, sg.z1_u], sg.is_os_sep * sg.dx) - sg.x2, sg.y2, sg.z2 = np.subtract([sg.x2_u, sg.y2_u, sg.z2_u], sg.is_os_sep * sg.dx) + # discretisted coordinates of the IS + sg.x1, sg.y1, sg.z1 = uip.round_to_grid(p1) + sg.x2, sg.y2, sg.z2 = uip.round_to_grid(p2) def set_working_region_cells(self, sg): """Number of cells in each dimension for the working region.""" @@ -77,6 +73,9 @@ class SubGridBase(UserObjectMulti): """Set number of iterations that will take place in the subgrid.""" sg.iterations = main.iterations * sg.ratio + def set_name(self, sg): + sg.name = self.kwargs['ID'] + def setup(self, sg, grid, uip): """"Common setup to both all subgrid types.""" p1 = self.kwargs['p1'] @@ -84,6 +83,8 @@ class SubGridBase(UserObjectMulti): p1, p2 = uip.check_box_points(p1, p2, self.__str__()) + self.set_name(sg) + self.check_filters(grid) self.set_discretisation(sg, grid) @@ -97,16 +98,6 @@ class SubGridBase(UserObjectMulti): # set the indices related to the subgrids main grid placement self.set_main_grid_indices(sg, grid, uip, p1, p2) - """ - try: - uip.check_box_points([sg.i0, sg.j0, sg.k0], - [sg.i1, sg.j1, sg.k1], cmd_str) - except CmdInputError: - es_f = 'The subgrid should extend at least {} cells' - es = es_f.format(sg.is_os_sep * 2) - raise CmdInputError(cmd_str, es) - """ - self.set_working_region_cells(sg) self.set_total_cells(sg) self.set_iterations(sg, grid) @@ -142,11 +133,9 @@ class SubGridHSG(SubGridBase): ratio=3, ID='', is_os_sep=3, - pml_separation=4, subgrid_pml_thickness=6, - interpolation='linear', + interpolation=3, loss_mechanism=False, - loss_factor=False, filter=True, **kwargs): """Constructor.""" @@ -157,12 +146,10 @@ class SubGridHSG(SubGridBase): kwargs['ratio'] = ratio kwargs['ID'] = ID kwargs['is_os_sep'] = is_os_sep - kwargs['pml_separation'] = pml_separation + kwargs['pml_separation'] = ratio // 2 + 2 kwargs['subgrid_pml_thickness'] = subgrid_pml_thickness kwargs['interpolation'] = interpolation kwargs['filter'] = filter - kwargs['loss_mechanism'] = loss_mechanism - kwargs['loss_factor'] = loss_factor super().__init__(**kwargs) self.order = 18 @@ -171,7 +158,7 @@ class SubGridHSG(SubGridBase): def create(self, grid, uip): sg = SubGridHSGUser(**self.kwargs) self.setup(sg, grid, uip) - if grid.messages: + if config.general['messages']: print(sg) return sg diff --git a/gprMax/updates.py b/gprMax/updates.py index a5e74388..b3ace83d 100644 --- a/gprMax/updates.py +++ b/gprMax/updates.py @@ -16,126 +16,134 @@ # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . from importlib import import_module +from gprMax.fields_outputs import store_outputs +import gprMax.config as config +from gprMax.cython.fields_updates_normal import update_electric +from gprMax.cython.fields_updates_normal import update_magnetic + class CPUUpdates: def __init__(self, G): - self.G = G + self.grid = G self.dispersive_update_a = None self.dispersive_update_b = None def store_outputs(self, iteration): # Store field component values for every receiver and transmission line store_outputs(iteration, - self.G.Ex, - self.G.Ey, - self.G.Ez, - self.G.Hx, - self.G.Hy, - self.G.Hz, - self.G) + self.grid.Ex, + self.grid.Ey, + self.grid.Ez, + self.grid.Hx, + self.grid.Hy, + self.grid.Hz, + self.grid) def store_snapshots(self, iteration): # Store any snapshots - for snap in self.G.snapshots: + for snap in self.grid.snapshots: if snap.time == iteration + 1: - snap.store(self.G) + snap.store(self.grid) def update_magnetic(self): # Update magnetic field components - update_magnetic(self.G.nx, - self.G.ny, - self.G.nz, + update_magnetic(self.grid.nx, + self.grid.ny, + self.grid.nz, config.hostinfo['ompthreads'], - self.G.updatecoeffsH, - self.G.ID, - self.G.Ex, - self.G.Ey, - self.G.Ez, - self.G.Hx, - self.G.Hy, - self.G.Hz) + self.grid.updatecoeffsH, + self.grid.ID, + self.grid.Ex, + self.grid.Ey, + self.grid.Ez, + self.grid.Hx, + self.grid.Hy, + self.grid.Hz) - def update_magnetic_pml(self, iteration): + def update_magnetic_pml(self): # Update magnetic field components with the PML correction - for pml in self.G.pmls: - pml.update_magnetic(self.G) + for pml in self.grid.pmls: + pml.update_magnetic(self.grid) def update_magnetic_sources(self, iteration): # Update magnetic field components from sources - for source in self.G.transmissionlines + self.G.magneticdipoles: + for source in self.grid.transmissionlines + self.grid.magneticdipoles: source.update_magnetic(iteration, - self.G.updatecoeffsH, - self.G.ID, - self.G.Hx, - self.G.Hy, - self.G.Hz, - self.G) + self.grid.updatecoeffsH, + self.grid.ID, + self.grid.Hx, + self.grid.Hy, + self.grid.Hz, + self.grid) def update_electric_a(self): # Update electric field components # All materials are non-dispersive so do standard update - if Material.maxpoles == 0: - update_electric(self.G.nx, - self.G.ny, - self.G.nz, + if config.materials['maxpoles'] == 0: + update_electric(self.grid.nx, + self.grid.ny, + self.grid.nz, config.hostinfo['ompthreads'], - self.G.updatecoeffsE, - self.G.ID, - self.G.Ex, - self.G.Ey, - self.G.Ez, - self.G.Hx, - self.G.Hy, - self.G.Hz) + self.grid.updatecoeffsE, + self.grid.ID, + self.grid.Ex, + self.grid.Ey, + self.grid.Ez, + self.grid.Hx, + self.grid.Hy, + self.grid.Hz) # If there are any dispersive materials do 1st part of dispersive update # (it is split into two parts as it requires present and updated electric field values). - self.dispersive_update_a(self.G.nx, - self.G.ny, - self.G.nz, - config.hostinfo['ompthreads'], - self.G.updatecoeffsE, - self.G.updatecoeffsdispersive, - self.G.ID, - self.G.Tx, - self.G.Ty, - self.G.Tz, - self.G.Ex, - self.G.Ey, - self.G.Ez, - self.G.Hx, - self.G.Hy, - self.G.Hz) + else: + self.dispersive_update_a(self.grid.nx, + self.grid.ny, + self.grid.nz, + config.hostinfo['ompthreads'], + config.materials['maxpoles'], + self.grid.updatecoeffsE, + self.grid.updatecoeffsdispersive, + self.grid.ID, + self.grid.Tx, + self.grid.Ty, + self.grid.Tz, + self.grid.Ex, + self.grid.Ey, + self.grid.Ez, + self.grid.Hx, + self.grid.Hy, + self.grid.Hz) def update_electric_pml(self): # Update electric field components with the PML correction - for pml in self.G.pmls: - pml.update_electric(self.G) + for pml in self.grid.pmls: + pml.update_electric(self.grid) def update_electric_sources(self, iteration): # Update electric field components from sources (update any Hertzian dipole sources last) - for source in self.G.voltagesources + self.G.transmissionlines + self.G.hertziandipoles: - source.update_electric(iteration, self.G.updatecoeffsE, self.G.ID, self.G.Ex, self.G.Ey, self.G.Ez, self.G) + for source in self.grid.voltagesources + self.grid.transmissionlines + self.grid.hertziandipoles: + source.update_electric(iteration, self.grid.updatecoeffsE, self.grid.ID, self.grid.Ex, self.grid.Ey, self.grid.Ez, self.grid) 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. - update_e_dispersive_b(self.G.nx, - self.G.ny, - self.G.nz, - config.hostinfo['ompthreads'], - Material.maxpoles, - self.G.updatecoeffsdispersive, - self.G.ID, - self.G.Tx, - self.G.Ty, - self.G.Tz, - self.G.Ex, - self.G.Ey, - self.G.Ez) + if config.materials['maxpoles'] != 0: + self.dispersive_update_b(self.grid.nx, + self.grid.ny, + self.grid.nz, + config.hostinfo['ompthreads'], + config.materials['maxpoles'], + self.grid.updatecoeffsdispersive, + self.grid.ID, + self.grid.Tx, + self.grid.Ty, + self.grid.Tz, + self.grid.Ex, + self.grid.Ey, + self.grid.Ez) def adapt_dispersive_config(self, config): @@ -180,8 +188,5 @@ class CPUUpdates: self.dispersive_update_b = disp_b_f -class SubgridUpdates: - pass - class GPUUpdates: pass diff --git a/gprMax/user_inputs.py b/gprMax/user_inputs.py index c7a23091..b269858e 100644 --- a/gprMax/user_inputs.py +++ b/gprMax/user_inputs.py @@ -1,6 +1,8 @@ from .exceptions import CmdInputError from .subgrids.base import SubGridBase from .utilities import round_value +from .utilities import Printer +import gprMax.config as config import numpy as np from colorama import init @@ -16,11 +18,15 @@ 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): """Return a point checker class based on the grid supplied.""" if isinstance(grid, SubGridBase): - return SubgridUserInput(grid) + + if config.general['autotranslate']: + print('AUTO-TRANSLATE SUB-GRID POINT ON') + return SubgridUserInput(grid) + else: + return MainGridUserInput(grid) else: return MainGridUserInput(grid) @@ -102,16 +108,18 @@ class MainGridUserInput(UserInput): class SubgridUserInput(MainGridUserInput): - """Class to handle (x, y, z) points supplied by the user in the sub grid.""" + """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.""" def __init__(self, grid): super().__init__(grid) # defines the region exposed to the user self.inner_bound = np.array([ - grid.pmlthickness['x0'] + grid.pml_separation, - grid.pmlthickness['y0'] + grid.pml_separation, - grid.pmlthickness['z0'] + grid.pml_separation]) + grid.n_boundary_cells_x, + grid.n_boundary_cells_y, + grid.n_boundary_cells_z]) self.outer_bound = np.subtract([grid.nx, grid.ny, grid.nz], self.inner_bound) @@ -119,12 +127,18 @@ class SubgridUserInput(MainGridUserInput): def translate_to_gap(self, p): """Function to translate the user input point to the real point in the subgrid""" - return np.add(p, self.inner_bound) + + 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 + p3 = (p[2] - self.grid.k0 * self.grid.ratio) + self.grid.n_boundary_cells_z + + 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""" + p = super().discretise_point(p) p_t = self.translate_to_gap(p) return p_t diff --git a/gprMax/utilities.py b/gprMax/utilities.py index d9a8cc30..ef92238e 100644 --- a/gprMax/utilities.py +++ b/gprMax/utilities.py @@ -418,3 +418,12 @@ def detect_check_gpus(deviceIDs): def timer(): """Function to return the current process wide time in fractional seconds.""" return perf_counter() + +class Printer(): + + def __init__(self, sim_config): + self.printing = sim_config.general['messages'] + + def print(self, str): + if self.printing: + print(str)