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)