subgridding files

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'])
+            weighting = np.array(self.kwargs['weighting'], dtype=np.float64)
-            limits = self.kwargs['limits']
+            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))

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:

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

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,
+#   autotranslate: auto translate objects with main grid coordinates
-           'progressbars': True, 'mode': '3D', 'cpu': True, 'cuda': False, 'opencl': False}
+#   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_output_file_path()
-        self.set_model_end()
+        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
+        # set range for number of models to run (internally 0 index)
-        if not self.mpi and not self.mpi_no_spawn:
+        if self.args.task:
-            # Set range for number of models to run
+            # Job array feeds args.n number of single tasks
-            if self.args.task:
+            modelstart = self.args.task - 1
-              # Job array feeds args.n number of single tasks
+            modelend = self.args.task
-              self.model_start = self.args.task - 1
+        elif self.args.restart:
-            elif self.args.restart:
+            modelstart = self.args.restart - 1
-              self.model_start = self.args.restart - 1
+            modelend = modelstart + self.args.n - 1
-            else:
+        else:
-              self.model_start = 0
+            modelstart = 0
+            modelend = modelstart + self.args.n
 
-        # mpi simulation
+        self.model_start = modelstart
-        elif self.mpi:
+        self.model_end = modelend
-            # 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!
 
     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

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

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:

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

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.

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)
-    context = create_context(sim_config, solver)
     context.run()

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

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

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 = SoilPeplinski(sand_fraction=float(tmp[0]),
-            soil = PeplinskiSoil(
+                                 clay_fraction=float(tmp[1]),
-                sand_fraction=float(tmp[0]),
+                                 bulk_density=float(tmp[2]),
-                clay_fraction=float(tmp[1]),
+                                 sand_density=float(tmp[3]),
-                bulk_density=float(tmp[2]),
+                                 water_fraction_lower=float(tmp[4]),
-                sand_density=float(tmp[3]),
+                                 water_fraction_upper=float(tmp[5]),
-                water_fraction_lower=float(tmp[4]),
+                                 ID=tmp[6])
-                water_fraction_upper=float(tmp[5]),
-                ID=tmp[6])
             scene_objects.append(soil)
 
     # Geometry views (creates VTK-based geometry files)

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):
+    def __init__(self, G, sim_config, model_config):
-        self.solver = solver
+        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:
+        if self.model_config.reuse_geometry:
-            self.build_geometry()
-        else:
             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
+        # combine available grids
-        printer.print(pml_information(G))
+        grids = [G] + G.subgrids
+        gridbuilders = [GridBuilder(grid, self.printer) for grid in grids]
 
-        self.build_pmls()
+        for gb in gridbuilders:
-        self.build_components()
+            gb.printer.print(pml_information(gb.grid))
-
+            gb.build_pmls()
-        # update grid for tm modes
+            gb.build_components()
-        self.tm_grid_update()
+            gb.tm_grid_update()
-
+            gb.update_voltage_source_materials()
-        self.update_voltage_source_materials()
+            gb.grid.initialise_std_update_coeff_arrays()
-
-        # Initialise arrays of update coefficients to pass to update functions
-        G.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,
+        for gb in gridbuilders:
-        # store in arrays, and build text list of materials/properties
+            gb.build_materials()
-        materialsdata = process_materials(G)
-        materialstable = SingleTable(materialsdata)
-        materialstable.outer_border = False
-        materialstable.justify_columns[0] = 'right'
-
-        printer.print('\nMaterials:')
-        printer.print(materialstable.table)
 
         # 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)
+        G = self.G
-        inputfilestr = '\n--- Model {}/{}, input file (not re-processed, i.e. geometry fixed): {}'.format(currentmodelrun, modelend, model_config.input_file_path)
+        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)
-        printer.print(Fore.GREEN + '{} {}\n'.format(model_config.inputfilestr, '-' * (get_terminal_width() - 1 - len(model_config.inputfilestr))) + Style.RESET_ALL)
+        self.printer.print(Fore.GREEN + '{} {}\n'.format(self.model_config.inputfilestr, '-' * (get_terminal_width() - 1 - len(inputfilestr))) + Style.RESET_ALL)
-        self.G.reset_fields()
+        for grid in [G] + G.subgrids:
-
+            grid.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()
 
     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):
+    def create_output_directory(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'])
 
-        for pml_id, thickness in self.G.pmlthickness.items():
+        if self.G.outputdirectory:
-            build_pml(self.G, pml_id, thickness)
+            # 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
+        self.printer.print('')
-        printer = Printer(self.sim_config)
+        pbar = tqdm(total=2, desc='Building {} grid'.format(self.grid.name), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars'])
-        printer.print('')
+        build_electric_components(self.grid.solid, self.grid.rigidE, self.grid.ID, self.grid)
-        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)
         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:
+        for voltagesource in self.grid.voltagesources:
-            voltagesource.create_material(self.G)
+            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):
+        self.printer.print('\n{} Grid Materials:'.format(self.grid.name))
-
+        self.printer.print(materialstable.table)
-        # 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

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.children_multiple, sg)
-            self.process_cmds(sg_cmd.geometry_cmds, sg)
+            self.process_cmds(sg_cmd.children_geometry, sg)
 
     def process_cmds(self, commands, grid):
         cmds_sorted = sorted(commands, key=lambda cmd: cmd.order)

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)
+    elif sim_config.subgrid:
-    else:
-        from .grid import FDTDGrid
         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()

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

gprMax/subgrids/base.py

@@ -17,8 +17,6 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 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.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_p = np.array([self.n_boundary_cells_x, self.n_boundary_cells_y, self.n_boundary_cells_z])
-                                           self.n_boundary_cells_y, self.n_boundary_cells_z)
 
         # interpolation scheme
         self.interpolation = kwargs['interpolation']

gprMax/subgrids/factory.py

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

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,

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

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

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.dl = grid.dl / sg.ratio
-        sg.dy = grid.dy / sg.ratio
+        sg.dx, sg.dy, sg.dz = grid.dl
-        sg.dz = grid.dz / sg.ratio
 
     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
+        # IS indices
-        sg.i0, sg.j0, sg.k0 = np.add([sg.i0_u, sg.j0_u, sg.k0_u], sg.is_os_sep)
+        sg.i0, sg.j0, sg.k0 = p1
-        sg.i1, sg.j1, sg.k1 = np.subtract([sg.i1_u, sg.j1_u, sg.k1_u], sg.is_os_sep)
+        sg.i1, sg.j1, sg.k1 = p2
 
-        # Main grid indices of the sub grid. These are dummy indices. They are
+        # OS indices
-        # not user internal except for printing to the user
+        sg.i_l, sg.j_l, sg.k_l = p1 - sg.is_os_sep
-        sg.x1_u, sg.y1_u, sg.z1_u = uip.round_to_grid(p1)
+        sg.i_u, sg.j_u, sg.k_u = p2 + sg.is_os_sep
-        sg.x2_u, sg.y2_u, sg.z2_u = uip.round_to_grid(p2)
 
-        sg.x1, sg.y1, sg.z1 = np.add([sg.x1_u, sg.y1_u, sg.z1_u], sg.is_os_sep * sg.dx)
+        # discretisted coordinates of the IS
-        sg.x2, sg.y2, sg.z2 = np.subtract([sg.x2_u, sg.y2_u, sg.z2_u], sg.is_os_sep * sg.dx)
+        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
 

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 <http://www.gnu.org/licenses/>.
 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.grid.Ex,
-                      self.G.Ey,
+                      self.grid.Ey,
-                      self.G.Ez,
+                      self.grid.Ez,
-                      self.G.Hx,
+                      self.grid.Hx,
-                      self.G.Hy,
+                      self.grid.Hy,
-                      self.G.Hz,
+                      self.grid.Hz,
-                      self.G)
+                      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,
+        update_magnetic(self.grid.nx,
-                        self.G.ny,
+                        self.grid.ny,
-                        self.G.nz,
+                        self.grid.nz,
                         config.hostinfo['ompthreads'],
-                        self.G.updatecoeffsH,
+                        self.grid.updatecoeffsH,
-                        self.G.ID,
+                        self.grid.ID,
-                        self.G.Ex,
+                        self.grid.Ex,
-                        self.G.Ey,
+                        self.grid.Ey,
-                        self.G.Ez,
+                        self.grid.Ez,
-                        self.G.Hx,
+                        self.grid.Hx,
-                        self.G.Hy,
+                        self.grid.Hy,
-                        self.G.Hz)
+                        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:
+        for pml in self.grid.pmls:
-            pml.update_magnetic(self.G)
+            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.grid.updatecoeffsH,
-                                   self.G.ID,
+                                   self.grid.ID,
-                                   self.G.Hx,
+                                   self.grid.Hx,
-                                   self.G.Hy,
+                                   self.grid.Hy,
-                                   self.G.Hz,
+                                   self.grid.Hz,
-                                   self.G)
+                                   self.grid)
 
     def update_electric_a(self):
         # Update electric field components
         # All materials are non-dispersive so do standard update
-        if Material.maxpoles == 0:
+        if config.materials['maxpoles'] == 0:
-            update_electric(self.G.nx,
+            update_electric(self.grid.nx,
-                            self.G.ny,
+                            self.grid.ny,
-                            self.G.nz,
+                            self.grid.nz,
                             config.hostinfo['ompthreads'],
-                            self.G.updatecoeffsE,
+                            self.grid.updatecoeffsE,
-                            self.G.ID,
+                            self.grid.ID,
-                            self.G.Ex,
+                            self.grid.Ex,
-                            self.G.Ey,
+                            self.grid.Ey,
-                            self.G.Ez,
+                            self.grid.Ez,
-                            self.G.Hx,
+                            self.grid.Hx,
-                            self.G.Hy,
+                            self.grid.Hy,
-                            self.G.Hz)
+                            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,
+        else:
-                                   self.G.ny,
+            self.dispersive_update_a(self.grid.nx,
-                                   self.G.nz,
+                                     self.grid.ny,
-                                   config.hostinfo['ompthreads'],
+                                     self.grid.nz,
-                                   self.G.updatecoeffsE,
+                                     config.hostinfo['ompthreads'],
-                                   self.G.updatecoeffsdispersive,
+                                     config.materials['maxpoles'],
-                                   self.G.ID,
+                                     self.grid.updatecoeffsE,
-                                   self.G.Tx,
+                                     self.grid.updatecoeffsdispersive,
-                                   self.G.Ty,
+                                     self.grid.ID,
-                                   self.G.Tz,
+                                     self.grid.Tx,
-                                   self.G.Ex,
+                                     self.grid.Ty,
-                                   self.G.Ey,
+                                     self.grid.Tz,
-                                   self.G.Ez,
+                                     self.grid.Ex,
-                                   self.G.Hx,
+                                     self.grid.Ey,
-                                   self.G.Hy,
+                                     self.grid.Ez,
-                                   self.G.Hz)
+                                     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:
+        for pml in self.grid.pmls:
-            pml.update_electric(self.G)
+            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:
+        for source in self.grid.voltagesources + self.grid.transmissionlines + self.grid.hertziandipoles:
-            source.update_electric(iteration, self.G.updatecoeffsE, self.G.ID, self.G.Ex, self.G.Ey, self.G.Ez, self.G)
+            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,
+        if config.materials['maxpoles'] != 0:
-                            self.G.ny,
+            self.dispersive_update_b(self.grid.nx,
-                            self.G.nz,
+                                     self.grid.ny,
-                            config.hostinfo['ompthreads'],
+                                     self.grid.nz,
-                            Material.maxpoles,
+                                     config.hostinfo['ompthreads'],
-                            self.G.updatecoeffsdispersive,
+                                     config.materials['maxpoles'],
-                            self.G.ID,
+                                     self.grid.updatecoeffsdispersive,
-                            self.G.Tx,
+                                     self.grid.ID,
-                            self.G.Ty,
+                                     self.grid.Tx,
-                            self.G.Tz,
+                                     self.grid.Ty,
-                            self.G.Ex,
+                                     self.grid.Tz,
-                            self.G.Ey,
+                                     self.grid.Ex,
-                            self.G.Ez)
+                                     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

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.n_boundary_cells_x,
-                                    grid.pmlthickness['y0'] + grid.pml_separation,
+                                    grid.n_boundary_cells_y,
-                                    grid.pmlthickness['z0'] + grid.pml_separation])
+                                    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

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)