diff --git a/gprMax/cmds_geometry/add_grass.py b/gprMax/cmds_geometry/add_grass.py
index efc052a8..9277bfe7 100644
--- a/gprMax/cmds_geometry/add_grass.py
+++ b/gprMax/cmds_geometry/add_grass.py
@@ -160,5 +160,5 @@ class AddGrass(UserObjectGeometry):
volume.fractalsurfaces.append(surface)
- if grid.messages:
+ if config.is_messages():
tqdm.write('{} blades of grass on surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with fractal dimension {:g}, fractal seeding {}, and range {:g}m to {:g}m, added to {}.'.format(n_blades, xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, surface.dimension, surface.seed, limits[0], limits[1], surface.operatingonID))
diff --git a/gprMax/cmds_geometry/add_surface_roughness.py b/gprMax/cmds_geometry/add_surface_roughness.py
index 227110b8..bf35549e 100644
--- a/gprMax/cmds_geometry/add_surface_roughness.py
+++ b/gprMax/cmds_geometry/add_surface_roughness.py
@@ -123,5 +123,5 @@ class AddSurfaceRoughness(UserObjectGeometry):
surface.generate_fractal_surface(grid)
volume.fractalsurfaces.append(surface)
- if config.general['messages']:
+ if config.is_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/cmds_geometry/add_surface_water.py b/gprMax/cmds_geometry/add_surface_water.py
index 75e80629..b954e2bd 100644
--- a/gprMax/cmds_geometry/add_surface_water.py
+++ b/gprMax/cmds_geometry/add_surface_water.py
@@ -114,5 +114,5 @@ class AddSurfaceWater(UserObjectGeometry):
if testwater:
raise CmdInputError(self.__str__() + ' requires the time step for the model to be less than the relaxation time required to model water.')
- if grid.messages:
+ if config.is_messages():
tqdm.write('Water on surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with depth {:g}m, added to {}.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, filldepth, surface.operatingonID))
diff --git a/gprMax/cmds_geometry/box.py b/gprMax/cmds_geometry/box.py
index b82774c0..9fd674ed 100644
--- a/gprMax/cmds_geometry/box.py
+++ b/gprMax/cmds_geometry/box.py
@@ -3,6 +3,7 @@ from .cmds_geometry import UserObjectGeometry
from ..exceptions import CmdInputError
from ..materials import Material
from ..cython.geometry_primitives import build_box
+import gprMax.config as config
from tqdm import tqdm
import numpy as np
@@ -84,7 +85,7 @@ class Box(UserObjectGeometry):
build_box(xs, xf, ys, yf, zs, zf, numID, numIDx, numIDy, numIDz, averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
- if grid.messages:
+ if config.is_messages():
if averaging:
dielectricsmoothing = 'on'
else:
diff --git a/gprMax/cmds_geometry/cylinder.py b/gprMax/cmds_geometry/cylinder.py
index c66ce2dc..ec1ea061 100644
--- a/gprMax/cmds_geometry/cylinder.py
+++ b/gprMax/cmds_geometry/cylinder.py
@@ -6,6 +6,7 @@ from ..cython.geometry_primitives import build_cylinder
from tqdm import tqdm
import numpy as np
+import gprMax.config as config
class Cylinder(UserObjectGeometry):
@@ -89,7 +90,7 @@ class Cylinder(UserObjectGeometry):
build_cylinder(x1, y1, z1, x2, y2, z2, r, grid.dx, grid.dy, grid.dz, numID, numIDx, numIDy, numIDz, averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
- if grid.messages:
+ if config.is_messages():
if averaging:
dielectricsmoothing = 'on'
else:
diff --git a/gprMax/cmds_geometry/cylindrical_sector.py b/gprMax/cmds_geometry/cylindrical_sector.py
index 3c71769d..c572417c 100644
--- a/gprMax/cmds_geometry/cylindrical_sector.py
+++ b/gprMax/cmds_geometry/cylindrical_sector.py
@@ -126,7 +126,7 @@ class CylindricalSector(UserObjectGeometry):
build_cylindrical_sector(ctr1, ctr2, level, sectorstartangle, sectorangle, r, normal, thickness, grid.dx, grid.dy, grid.dz, numID, numIDx, numIDy, numIDz, averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
- if grid.messages:
+ if config.is_messages():
if thickness > 0:
if averaging:
dielectricsmoothing = 'on'
diff --git a/gprMax/cmds_geometry/edge.py b/gprMax/cmds_geometry/edge.py
index fd01b97a..1e4ff282 100644
--- a/gprMax/cmds_geometry/edge.py
+++ b/gprMax/cmds_geometry/edge.py
@@ -60,5 +60,5 @@ class Edge(UserObjectGeometry):
for k in range(zs, zf):
build_edge_z(xs, ys, k, material.numID, grid.rigidE, grid.rigidH, grid.ID)
- if grid.messages:
+ if config.is_messages():
tqdm.write('Edge from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material {} created.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, material_id))
diff --git a/gprMax/cmds_geometry/fractal_box.py b/gprMax/cmds_geometry/fractal_box.py
index 959a37dd..50ec3b63 100644
--- a/gprMax/cmds_geometry/fractal_box.py
+++ b/gprMax/cmds_geometry/fractal_box.py
@@ -25,7 +25,7 @@ class FractalBox(UserObjectGeometry):
weighting = np.array(self.kwargs['weighting'])
n_materials = self.kwargs['n_materials']
mixing_model_id = self.kwargs['mixing_model_id']
- ID = self.kwargs['ID']
+ ID = self.kwargs['id']
except KeyError:
raise CmdInputError(self.__str__() + ' Incorrect parameters')
@@ -83,7 +83,7 @@ class FractalBox(UserObjectGeometry):
volume.averaging = averagefractalbox
volume.mixingmodel = mixingmodel
- if config.general['messages']:
+ if config.is_messages():
if volume.averaging:
dielectricsmoothing = 'on'
else:
diff --git a/gprMax/cmds_geometry/plate.py b/gprMax/cmds_geometry/plate.py
index a2b63340..74ebaed0 100644
--- a/gprMax/cmds_geometry/plate.py
+++ b/gprMax/cmds_geometry/plate.py
@@ -108,5 +108,5 @@ class Plate(UserObjectGeometry):
for j in range(ys, yf):
build_face_xy(i, j, zs, numIDx, numIDy, grid.rigidE, grid.rigidH, grid.ID)
- if grid.messages:
+ if config.is_messages():
tqdm.write('Plate from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material(s) {} created.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, ', '.join(materialsrequested)))
diff --git a/gprMax/cmds_geometry/sphere.py b/gprMax/cmds_geometry/sphere.py
index 590a7ad9..e38ea409 100644
--- a/gprMax/cmds_geometry/sphere.py
+++ b/gprMax/cmds_geometry/sphere.py
@@ -83,7 +83,7 @@ class Sphere(UserObjectGeometry):
build_sphere(xc, yc, zc, r, grid.dx, grid.dy, grid.dz, numID, numIDx, numIDy, numIDz, averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
- if grid.messages:
+ if config.is_messages():
if averaging:
dielectricsmoothing = 'on'
else:
diff --git a/gprMax/cmds_geometry/triangle.py b/gprMax/cmds_geometry/triangle.py
index bc68c9c2..be2eadd9 100644
--- a/gprMax/cmds_geometry/triangle.py
+++ b/gprMax/cmds_geometry/triangle.py
@@ -121,7 +121,7 @@ class Triangle(UserObjectGeometry):
build_triangle(x1, y1, z1, x2, y2, z2, x3, y3, z3, normal, thickness, grid.dx, grid.dy, grid.dz, numID, numIDx, numIDy, numIDz, averaging, grid.solid, grid.rigidE, grid.rigidH, grid.ID)
- if grid.messages:
+ if config.is_messages():
if thickness > 0:
if averaging:
dielectricsmoothing = 'on'
diff --git a/gprMax/cmds_multiple.py b/gprMax/cmds_multiple.py
index 28a16dd1..8316c143 100644
--- a/gprMax/cmds_multiple.py
+++ b/gprMax/cmds_multiple.py
@@ -106,7 +106,7 @@ class Waveform(UserObjectMulti):
w.amp = amp
w.freq = freq
- if config.general['messages']:
+ if config.is_messages():
print('Waveform {} of type {} with maximum amplitude scaling {:g}, frequency {:g}Hz created.'.format(w.ID, w.type, w.amp, w.freq))
grid.waveforms.append(w)
@@ -183,7 +183,7 @@ class VoltageSource(UserObjectMulti):
v.calculate_waveform_values(grid)
- if config.general['messages']:
+ if config.is_messages():
print('Voltage source with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(v.polarisation, v.xcoord * grid.dx, v.ycoord * grid.dy, v.zcoord * grid.dz, v.resistance) + startstop + 'using waveform {} created.'.format(v.waveformID))
grid.voltagesources.append(v)
@@ -268,7 +268,7 @@ class HertzianDipole(UserObjectMulti):
h.calculate_waveform_values(grid)
- if config.general['messages']:
+ if config.is_messages():
if grid.mode == '2D':
print('Hertzian dipole is a line source in 2D with polarity {} at {:g}m, {:g}m, {:g}m,'.format(h.polarisation, h.xcoord * grid.dx, h.ycoord * grid.dy, h.zcoord * grid.dz) + startstop + 'using waveform {} created.'.format(h.waveformID))
else:
@@ -345,7 +345,7 @@ class MagneticDipole(UserObjectMulti):
m.calculate_waveform_values(grid)
- if config.general['messages']:
+ if config.is_messages():
print('Magnetic dipole with polarity {} at {:g}m, {:g}m, {:g}m,'.format(m.polarisation, m.xcoord * grid.dx, m.ycoord * grid.dy, m.zcoord * grid.dz) + startstop + 'using waveform {} created.'.format(m.waveformID))
grid.magneticdipoles.append(m)
@@ -427,7 +427,7 @@ class TransmissionLine(UserObjectMulti):
t.calculate_waveform_values(grid)
t.calculate_incident_V_I(grid)
- if config.general['messages']:
+ if config.is_messages():
print('Transmission line with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(t.polarisation, t.xcoord * grid.dx, t.ycoord * grid.dy, t.zcoord * grid.dz, t.resistance) + startstop + 'using waveform {} created.'.format(t.waveformID))
grid.transmissionlines.append(t)
@@ -456,7 +456,7 @@ class Rx(UserObjectMulti):
r.xcoordorigin, r.ycoordorigin, r.zcoordorigin = p
try:
- r.ID = self.kwargs['ID']
+ r.ID = self.kwargs['id']
outputs = self.kwargs['outputs']
# Get allowable outputs
if grid.gpu is not None:
@@ -475,8 +475,7 @@ class Rx(UserObjectMulti):
r.ID = r.__class__.__name__ + '(' + str(r.xcoord) + ',' + str(r.ycoord) + ',' + str(r.zcoord) + ')'
for key in RxUser.defaultoutputs:
r.outputs[key] = np.zeros(grid.iterations, dtype=floattype)
-
- if config.general['messages']:
+ if config.is_messages():
print('Receiver at {:g}m, {:g}m, {:g}m with output component(s) {} created.'.format(r.xcoord * grid.dx, r.ycoord * grid.dy, r.zcoord * grid.dz, ', '.join(r.outputs)))
grid.rxs.append(r)
@@ -526,7 +525,7 @@ class RxArray(UserObjectMulti):
else:
raise CmdInputError("'{}' the step size should not be less than the spatial discretisation".format(self.__str__()))
- if config.general['messages']:
+ if config.is_messages():
print('Receiver array {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with steps {:g}m, {:g}m, {:g}m'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, dx * grid.dx, dy * grid.dy, dz * grid.dz))
for x in range(xs, xf + 1, dx):
@@ -542,7 +541,7 @@ class RxArray(UserObjectMulti):
r.ID = r.__class__.__name__ + '(' + str(x) + ',' + str(y) + ',' + str(z) + ')'
for key in RxUser.defaultoutputs:
r.outputs[key] = np.zeros(grid.iterations, dtype=floattype)
- if config.general['messages']:
+ if config.is_messages():
print(' Receiver at {:g}m, {:g}m, {:g}m with output component(s) {} created.'.format(r.xcoord * grid.dx, r.ycoord * grid.dy, r.zcoord * grid.dz, ', '.join(r.outputs)))
grid.rxs.append(r)
@@ -601,7 +600,7 @@ class Snapshot(UserObjectMulti):
else:
s = SnapshotUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, iterations, filename)
- if config.general['messages']:
+ if config.is_messages():
print('Snapshot from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, at {:g} secs with filename {} created.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, dx * grid.dx, dy * grid.dy, dz * grid.dz, s.time * grid.dt, s.basefilename))
grid.snapshots.append(s)
@@ -651,7 +650,7 @@ class Material(UserObjectMulti):
if m.se == float('inf'):
m.averagable = False
- if config.general['messages']:
+ if config.is_messages():
tqdm.write('Material {} with eps_r={:g}, sigma={:g} S/m; mu_r={:g}, sigma*={:g} Ohm/m created.'.format(m.ID, m.er, m.se, m.mr, m.sm))
# Append the new material object to the materials list
@@ -702,7 +701,7 @@ class AddDebyeDispersion(UserObjectMulti):
if material.poles > MaterialUser.maxpoles:
MaterialUser.maxpoles = material.poles
- if config.general['messages']:
+ if config.is_messages():
tqdm.write('Debye disperion added to {} with delta_eps_r={}, and tau={} secs created.'.format(material.ID, ', '.join('%4.2f' % deltaer for deltaer in material.deltaer), ', '.join('%4.3e' % tau for tau in material.tau)))
@@ -750,7 +749,7 @@ class AddLorentzDispersion(UserObjectMulti):
if material.poles > MaterialUser.maxpoles:
MaterialUser.maxpoles = material.poles
- if config.general['messages']:
+ if config.is_messages():
tqdm.write('Lorentz disperion added to {} with delta_eps_r={}, omega={} secs, and gamma={} created.'.format(material.ID, ', '.join('%4.2f' % deltaer for deltaer in material.deltaer), ', '.join('%4.3e' % tau for tau in material.tau), ', '.join('%4.3e' % alpha for alpha in material.alpha)))
@@ -795,7 +794,7 @@ class AddDrudeDispersion(UserObjectMulti):
if material.poles > MaterialUser.maxpoles:
MaterialUser.maxpoles = material.poles
- if config.general['messages']:
+ if config.is_messages():
tqdm.write('Drude disperion added to {} with omega={} secs, and gamma={} secs created.'.format(material.ID, ', '.join('%4.3e' % tau for tau in material.tau), ', '.join('%4.3e' % alpha for alpha in material.alpha)))
@@ -817,7 +816,7 @@ class SoilPeplinski(UserObjectMulti):
sand_density = self.kwargs['sand_density']
water_fraction_lower = self.kwargs['water_fraction_lower']
water_fraction_upper = self.kwargs['water_fraction_upper']
- ID = self.kwargs['ID']
+ ID = self.kwargs['id']
except KeyError:
raise CmdInputError(self.__str__() + ' requires at exactly seven parameters')
@@ -840,7 +839,7 @@ class SoilPeplinski(UserObjectMulti):
# Create a new instance of the Material class material (start index after pec & free_space)
s = PeplinskiSoilUser(ID, sand_fraction, clay_fraction, bulk_density, sand_density, (water_fraction_lower, water_fraction_upper))
- if config.general['messages']:
+ if config.is_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]))
# Append the new material object to the materials list
@@ -912,7 +911,7 @@ class GeometryView(UserObjectMulti):
g = GeometryViewUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, fileext, grid)
- if config.general['messages']:
+ if config.is_messages():
print('Geometry view from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, multi_grid {}, grid={}, with filename base {} created.'.format(xs * grid.dx, ys * grid.dy, zs * grid.dz, xf * grid.dx, yf * grid.dy, zf * grid.dz, dx * grid.dx, dy * grid.dy, dz * grid.dz, self.multi_grid, grid.name, g.basefilename))
# Append the new GeometryView object to the geometry views list
@@ -942,7 +941,7 @@ class GeometryObjectsWrite(UserObjectMulti):
g = GeometryObjectsUser(x0, y0, z0, x1, y1, z1, filename)
- if config.general['messages']:
+ if config.is_messages():
print('Geometry objects in the volume from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, will be written to {}, with materials written to {}'.format(p1[0] * grid.dx, p1[1] * grid.dy, p1[2] * grid.dz, p2[0] * grid.dx, p2[1] * grid.dy, p2[2] * grid.dz, g.filename, g.materialsfilename))
# Append the new GeometryView object to the geometry objects to write list
@@ -1017,7 +1016,7 @@ class PMLCFS(UserObjectMulti):
cfs.kappa = cfskappa
cfs.sigma = cfssigma
- if config.general['messages']:
+ if config.is_messages():
print('PML CFS parameters: alpha (scaling: {}, scaling direction: {}, min: {:g}, max: {:g}), kappa (scaling: {}, scaling direction: {}, min: {:g}, max: {:g}), sigma (scaling: {}, scaling direction: {}, min: {:g}, max: {}) created.'.format(cfsalpha.scalingprofile, cfsalpha.scalingdirection, cfsalpha.min, cfsalpha.max, cfskappa.scalingprofile, cfskappa.scalingdirection, cfskappa.min, cfskappa.max, cfssigma.scalingprofile, cfssigma.scalingdirection, cfssigma.min, cfssigma.max))
grid.cfs.append(cfs)
diff --git a/gprMax/cmds_single_use.py b/gprMax/cmds_single_use.py
index 61d95eb4..fe1aa759 100644
--- a/gprMax/cmds_single_use.py
+++ b/gprMax/cmds_single_use.py
@@ -86,7 +86,7 @@ class DomainSingle(UserObjectSingle):
if G.nx == 0 or G.ny == 0 or G.nz == 0:
raise CmdInputError(self.__str__ + ' requires at least one cell in every dimension')
- if config.general['messages']:
+ if config.is_messages():
print('Domain size: {:g} x {:g} x {:g}m ({:d} x {:d} x {:d} = {:g} cells)'.format(self.kwargs['p1'][0], self.kwargs['p1'][1], self.kwargs['p1'][2], G.nx, G.ny, G.nz, (G.nx * G.ny * G.nz)))
# Time step CFL limit (either 2D or 3D); switch off appropriate PMLs for 2D
@@ -114,7 +114,7 @@ class DomainSingle(UserObjectSingle):
# Avoids inadvertently exceeding the CFL due to binary representation of floating point number.
G.dt = round_value(G.dt, decimalplaces=d.getcontext().prec - 1)
- if config.general['messages']:
+ if config.is_messages():
print('Mode: {}'.format(G.mode))
print('Time step (at CFL limit): {:g} secs'.format(G.dt))
@@ -139,7 +139,7 @@ class DomainSingle(UserObjectSingle):
G.nthreads = hostinfo['physicalcores']
os.environ['OMP_NUM_THREADS'] = str(G.nthreads)
- if config.general['messages']:
+ if config.is_messages():
print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
if G.nthreads > hostinfo['physicalcores']:
print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, hostinfo['physicalcores']) + Style.RESET_ALL)
@@ -194,7 +194,7 @@ class Discretisation(UserObjectSingle):
if G.dl[2] <= 0:
raise CmdInputError('Discretisation requires the z-direction spatial step to be greater than zero')
- if config.general['messages']:
+ if config.is_messages():
print('Spatial discretisation: {:g} x {:g} x {:g}m'.format(*G.dl))
@@ -243,7 +243,7 @@ class TimeWindow(UserObjectSingle):
if not G.timewindow:
raise CmdInputError('TimeWindow: Specify a time or number of iterations')
- if config.general['messages']:
+ if config.is_messages():
print('Time window: {:g} secs ({} iterations)'.format(G.timewindow, G.iterations))
@@ -291,7 +291,7 @@ class Title(UserObjectSingle):
except KeyError:
pass
- if config.general['messages']:
+ if config.is_messages():
print('Model title: {}'.format(G.title))
class NumThreads(UserObjectSingle):
@@ -321,13 +321,13 @@ class NumThreads(UserObjectSingle):
G.nthreads = n
os.environ['OMP_NUM_THREADS'] = str(G.nthreads)
- if config.general['messages']:
+ if config.is_messages():
print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
if G.nthreads > hostinfo['physicalcores']:
print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, hostinfo['physicalcores']) + Style.RESET_ALL)
# Print information about any GPU in use
- if config.general['messages']:
+ if config.is_messages():
if G.gpu is not None:
print('GPU solving using: {} - {}'.format(G.gpu.deviceID, G.gpu.name))
@@ -358,7 +358,7 @@ class TimeStepStabilityFactor(UserObjectSingle):
if f <= 0 or f > 1:
raise CmdInputError(self.__str__() + ' requires the value of the time step stability factor to be between zero and one')
G.dt = G.dt * f
- if config.general['messages']:
+ if config.is_messages():
print('Time step (modified): {:g} secs'.format(G.dt))
@@ -412,7 +412,7 @@ class SrcSteps(UserObjectSingle):
except KeyError:
raise CmdInputError('#src_steps: requires exactly three parameters')
# src_steps
- if config.general['messages']:
+ if config.is_messages():
print('Simple sources will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.srcsteps[0] * G.dx, G.srcsteps[1] * G.dy, G.srcsteps[2] * G.dz))
@@ -429,7 +429,7 @@ class RxSteps(UserObjectSingle):
G.rxsteps = uip.discretise_point(self.kwargs['p1'])
except KeyError:
raise CmdInputError('#rx_steps: requires exactly three parameters')
- if config.general['messages']:
+ if config.is_messages():
print('All receivers will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.rxsteps[0] * G.dx, G.rxsteps[1] * G.dy, G.rxsteps[2] * G.dz))
@@ -455,7 +455,7 @@ class ExcitationFile(UserObjectSingle):
if not os.path.isfile(excitationfile):
excitationfile = os.path.abspath(os.path.join(G.inputdirectory, excitationfile))
- if config.general['messages']:
+ if config.is_messages():
print('\nExcitation file: {}'.format(excitationfile))
# Get waveform names
@@ -495,7 +495,7 @@ class ExcitationFile(UserObjectSingle):
# Interpolate waveform values
w.userfunc = interpolate.interp1d(waveformtime, singlewaveformvalues, **kwargs)
- if config.general['messages']:
+ if config.is_messages():
print('User waveform {} created using {} and, if required, interpolation parameters (kind: {}, fill value: {}).'.format(w.ID, timestr, kwargs['kind'], kwargs['fill_value']))
G.waveforms.append(w)
diff --git a/gprMax/config.py b/gprMax/config.py
index 618a6da0..1e6879bf 100644
--- a/gprMax/config.py
+++ b/gprMax/config.py
@@ -51,6 +51,12 @@ z0 = np.sqrt(m0 / e0)
general = {'inputfilepath': 'gprMax', 'outputfilepath': 'gprMax', 'messages': True,
'progressbars': True, 'mode': '3D', 'cpu': True, 'cuda': False, 'opencl': False, 'autotranslate': False}
+
+def is_messages():
+ """Function to return messages."""
+ return general['messages']
+
+
# Store information about host machine
hostinfo = get_host_info()
diff --git a/gprMax/cython/fields_updates_hsg.pyx b/gprMax/cython/fields_updates_hsg.pyx
index 317ae303..7c56b42a 100644
--- a/gprMax/cython/fields_updates_hsg.pyx
+++ b/gprMax/cython/fields_updates_hsg.pyx
@@ -1,8 +1,24 @@
+# Copyright (C) 2015-2017: 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 .
import numpy as np
cimport numpy as np
from cython.parallel import prange
-
cpdef void cython_update_electric_os(
np.float64_t[:, :] updatecoeffsE,
np.uint32_t[:, :, :, :] ID,
@@ -208,7 +224,7 @@ cpdef void cython_update_is(
int nwx,
int nwy,
int nwz,
- double[:, :] updatecoeffsE,
+ np.float64_t[:, :] updatecoeffsE,
np.uint32_t[:, :, :, :] ID,
int n,
int offset,
@@ -216,9 +232,9 @@ cpdef void cython_update_is(
int nwm,
int nwn,
int face,
- double[:, :, :] field,
- double[:, :] inc_field_l,
- double[:, :] inc_field_u,
+ np.float64_t[:, :, :] field,
+ np.float64_t[:, :] inc_field_l,
+ np.float64_t[:, :] inc_field_u,
Py_ssize_t lookup_id,
int sign_l,
int sign_u,
diff --git a/gprMax/fields_outputs.py b/gprMax/fields_outputs.py
index 7a9bdd39..d0a61b33 100644
--- a/gprMax/fields_outputs.py
+++ b/gprMax/fields_outputs.py
@@ -17,6 +17,7 @@
# along with gprMax. If not, see .
from string import Template
+from pathlib import Path
import h5py
@@ -93,6 +94,53 @@ __global__ void store_outputs(int NRX, int iteration, const int* __restrict__ rx
def write_hdf5_outputfile(outputfile, G):
+ write_hdf5_main_grid_outputfile(outputfile, G)
+ write_hdf5_sub_grid_outputfile(outputfile, G)
+
+
+def write_hdf5_main_grid_outputfile(outputfile, G):
+ """Write an output file in HDF5 format.
+
+ Args:
+ outputfile (str): Name of the output file.
+ Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic field values.
+ G (class): Grid class instance - holds essential parameters describing the model.
+ """
+
+ write_data(outputfile, G)
+
+
+def write_hdf5_sub_grid_outputfile(outputfile, G):
+ """Write an output file in HDF5 format.
+
+ Args:
+ outputfile (str): Name of the output file.
+ Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic field values.
+ G (class): Grid class instance - holds essential parameters describing the model.
+ """
+
+ stem = outputfile.stem
+ suffix = outputfile.suffix
+ parent = outputfile.parent
+
+ for sg in G.subgrids:
+
+ # create an outputfile for each subgrid
+ fp = stem + '_' + sg.name + suffix
+ fp = parent / Path(fp)
+
+ f = write_data(fp, sg)
+
+ # write some additional meta data about the subgrid
+ f.attrs['is_os_sep'] = sg.is_os_sep
+ f.attrs['pml_separation'] = sg.pml_separation
+ f.attrs['subgrid_pml_thickness'] = sg.pmlthickness['x0']
+ f.attrs['filter'] = sg.filter
+ f.attrs['ratio'] = sg.ratio
+ f.attrs['interpolation'] = sg.interpolation
+
+
+def write_data(outputfile, G):
"""Write an output file in HDF5 format.
Args:
@@ -146,3 +194,5 @@ def write_hdf5_outputfile(outputfile, G):
for output in rx.outputs:
f['/rxs/rx' + str(rxindex + 1) + '/' + output] = rx.outputs[output]
+
+ return f
diff --git a/gprMax/input_cmds_geometry.py b/gprMax/input_cmds_geometry.py
index 8010f223..a25f0ec0 100644
--- a/gprMax/input_cmds_geometry.py
+++ b/gprMax/input_cmds_geometry.py
@@ -246,13 +246,13 @@ def process_geometrycmds(geometry):
ID = tmp[13]
# without seed
if len(tmp) == 14:
- fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, ID=ID, n_materials=n_materials)
+ fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, id=ID, n_materials=n_materials)
# with seed
elif len(tmp) == 15:
- fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, ID=ID, n_materials=n_materials, seed=tmp[14])
+ fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, id=ID, n_materials=n_materials, seed=tmp[14])
# user specified averaging
elif len(tmp) == 16:
- fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, ID=ID, n_materials=n_materials, seed=tmp[14], averaging=tmp[15].lower())
+ fb = FractalBox(p1=p1, p2=p2, frac_dim=frac_dim, weighting=weighting, mixing_model_id=mixing_model_id, id=ID, n_materials=n_materials, seed=tmp[14], averaging=tmp[15].lower())
else:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py
index 23aaac9d..5ead64de 100644
--- a/gprMax/input_cmds_multiuse.py
+++ b/gprMax/input_cmds_multiuse.py
@@ -143,11 +143,10 @@ def process_multicmds(multicmds):
tmp = cmdinstance.split()
if len(tmp) != 3 and len(tmp) < 5:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' has an incorrect number of parameters')
-
if len(tmp) == 3:
rx = Rx(p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])))
else:
- rx = Rx(p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])), ID=tmp[3], outputs=tmp[4:])
+ rx = Rx(p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])), id=tmp[3], outputs=tmp[4:])
scene_objects.append(rx)
# Receiver array
@@ -196,7 +195,7 @@ def process_multicmds(multicmds):
if len(tmp) != 5:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires exactly five parameters')
- material = Material(er=float(tmp[0]), se=float(tmp[1]), mr=float(tmp[2]), sm=float(tmp[3]), ID=tmp[4])
+ material = Material(er=float(tmp[0]), se=float(tmp[1]), mr=float(tmp[2]), sm=float(tmp[3]), id=tmp[4])
scene_objects.append(material)
cmdname = '#add_dispersion_debye'
@@ -274,7 +273,7 @@ def process_multicmds(multicmds):
sand_density=float(tmp[3]),
water_fraction_lower=float(tmp[4]),
water_fraction_upper=float(tmp[5]),
- ID=tmp[6])
+ 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 41f205e9..e9af5b99 100644
--- a/gprMax/model_build_run.py
+++ b/gprMax/model_build_run.py
@@ -81,7 +81,7 @@ class ModelBuildRun:
self.G = G
self.sim_config = sim_config
self.model_config = model_config
- self.printer = Printer(sim_config)
+ self.printer = Printer(config)
# Monitor memory usage
self.p = None
@@ -130,9 +130,9 @@ class ModelBuildRun:
receiver.zcoord = receiver.zcoordorigin + (currentmodelrun - 1) * G.rxsteps[2]
# Write files for any geometry views and geometry object outputs
- if not (G.geometryviews or G.geometryobjectswrite) and self.sim_config.geometry_only and config.general['messages']:
+ if not (G.geometryviews or G.geometryobjectswrite) and self.sim_config.geometry_only and config.is_messages():
print(Fore.RED + '\nWARNING: No geometry views or geometry objects to output found.' + Style.RESET_ALL)
- if config.general['messages']: print()
+ if config.is_messages(): print()
for i, geometryview in enumerate(G.geometryviews):
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'])
@@ -153,7 +153,7 @@ class ModelBuildRun:
G = self.G
- printer = Printer(sim_config)
+ printer = Printer(config)
printer.print(model_config.next_model)
scene = self.build_scene()
@@ -300,7 +300,7 @@ class ModelBuildRun:
"""
G = self.G
- if config.general['messages']:
+ if config.is_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:
diff --git a/gprMax/solvers.py b/gprMax/solvers.py
index 2f7d9024..26de9bb8 100644
--- a/gprMax/solvers.py
+++ b/gprMax/solvers.py
@@ -21,7 +21,7 @@ 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
+from .subgrids.solver import create_solver as create_subgrid_solver
def create_G(sim_config):
@@ -40,20 +40,20 @@ def create_solver(G, sim_config):
"""Returns the configured solver."""
if sim_config.gpu:
updates = GPUUpdates(G)
+ solver = Solver(updates)
elif sim_config.subgrid:
- updates = create_subgrid_updates(G)
+ solver = create_subgrid_solver(G)
else:
updates = CPUUpdates(G)
+ solver = Solver(updates)
- solver = Solver(updates)
+ return solver
# 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
+ #props = updates.adapt_dispersive_config(config)
+ #updates.set_dispersive_updates(props)
class Solver:
@@ -76,7 +76,6 @@ class Solver:
def solve(self, iterator):
"""Time step the FDTD model."""
tsolvestart = timer()
-
for iteration in iterator:
self.updates.store_outputs(iteration)
self.updates.store_snapshots(iteration)
diff --git a/gprMax/sources.py b/gprMax/sources.py
index 4ce3c26b..a1a9039f 100644
--- a/gprMax/sources.py
+++ b/gprMax/sources.py
@@ -169,7 +169,6 @@ class HertzianDipole(Source):
j = self.ycoord
k = self.zcoord
componentID = 'E' + self.polarisation
-
if self.polarisation == 'x':
Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * self.waveformvaluesJ[iteration] * self.dl * (1 / (G.dx * G.dy * G.dz))
diff --git a/gprMax/subgrids/base.py b/gprMax/subgrids/base.py
index 4bb0298e..7b0fb741 100644
--- a/gprMax/subgrids/base.py
+++ b/gprMax/subgrids/base.py
@@ -1,24 +1,18 @@
-# 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 ..grid import FDTDGrid
+
+from ..materials import Material
+
+from scipy.constants import mu_0
+from scipy.constants import epsilon_0
+from scipy.constants import c
+
import numpy as np
+from colorama import init
+from colorama import Fore
+from colorama import Style
+init()
+
class SubGridBase(FDTDGrid):
@@ -26,8 +20,6 @@ class SubGridBase(FDTDGrid):
super().__init__()
self.mode = '3D'
-
- # subgridding ratio
self.ratio = kwargs['ratio']
if self.ratio % 2 == 0:
@@ -60,17 +52,21 @@ class SubGridBase(FDTDGrid):
self.n_boundary_cells_y = d_to_pml + self.pmlthickness['y0']
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])
-
- # interpolation scheme
self.interpolation = kwargs['interpolation']
- def main_grid_index_to_subgrid_index(self, p):
- """
- Return the equivalent spatial index of the global (i, j, k) in the subgrid.
- Args:
- p (numpy.array): i, j, k indices of a point in the main grid.
- """
- return self.n_boundary_cells_p + (p - self.p0) * self.ratio
+ 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 main_grid_index_to_subgrid_index(self, i, j, k):
+ i_s = self.n_boundary_cells_x + (i - self.i0) * self.ratio
+ j_s = self.n_boundary_cells_y + (j - self.j0) * self.ratio
+ k_s = self.n_boundary_cells_z + (k - self.k0) * self.ratio
+
+ return (i_s, j_s, k_s)
+
+ def initialise_geometry_arrays(self):
+ super().initialise_geometry_arrays()
diff --git a/gprMax/subgrids/precursor_nodes.py b/gprMax/subgrids/precursor_nodes.py
index a57686ce..a837d18a 100644
--- a/gprMax/subgrids/precursor_nodes.py
+++ b/gprMax/subgrids/precursor_nodes.py
@@ -1,23 +1,6 @@
-# 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 .
-
import numpy as np
from scipy import interpolate
+import sys
def calculate_weighting_coefficients(x1, x):
@@ -26,14 +9,13 @@ def calculate_weighting_coefficients(x1, x):
return (c1, c2)
-class PrecusorNodesBase:
+class PrecusorNodes2dBase(object):
def __init__(self, fdtd_grid, sub_grid):
self.G = fdtd_grid
self.ratio = sub_grid.ratio
self.nwx = sub_grid.nwx
self.nwy = sub_grid.nwy
- self.nwz = sub_grid.nwz
self.sub_grid = sub_grid
self.interpolation = sub_grid.interpolation
@@ -60,6 +42,300 @@ class PrecusorNodesBase:
self.l_weight = self.ratio // 2
self.r_weight = self.ratio - self.l_weight
+ #self.l_weight = 1
+ #self.r_weight = 2
+
+ def _initialize_fields(self):
+ print('dont call me')
+ sys.exit()
+ pass
+
+ def _initialize_field_names(self):
+ print('dont call me')
+ sys.exit()
+ pass
+
+ def interpolate_magnetic_in_time(self, m):
+ self.weight_pre_and_current_fields(m, self.fn_m)
+
+ def interpolate_electric_in_time(self, m):
+ self.weight_pre_and_current_fields(m, self.fn_e)
+
+ def weight_pre_and_current_fields(self, m, field_names):
+ c1, c2 = calculate_weighting_coefficients(m, self.ratio)
+
+ for f in field_names:
+ try:
+ val = c1 * getattr(self, f + '_0') + c2 * getattr(self, f + '_1')
+ except ValueError:
+ print(self.ex_front_0.shape)
+ print(self.ex_front_1.shape)
+ raise Exception(f)
+ setattr(self, f, val)
+
+ def calc_exact_field(self, field_names):
+ """
+ Function to set the fields used in update calculations to the
+ values at the current main time step.
+ i.e. ey_left = copy.ey_left_1
+ """
+ for f in field_names:
+ val = np.copy(getattr(self, f + '_1'))
+ setattr(self, f, val)
+
+ def calc_exact_magnetic_in_time(self):
+ self.calc_exact_field(self.fn_m)
+
+ def calc_exact_electric_in_time(self):
+ self.calc_exact_field(self.fn_e)
+
+
+class PrecursorNodes2dTM(PrecusorNodes2dBase):
+
+ def __init__(self, fdtd_grid, sub_grid):
+ super().__init__(fdtd_grid, sub_grid)
+
+ def _initialize_fields(self):
+
+ # H precursors
+ self.hx_front_0 = np.zeros(self.nwx + 1)
+ self.hx_front_1 = np.zeros(self.nwx + 1)
+
+ self.hx_back_0 = np.zeros(self.nwx + 1)
+ self.hx_back_1 = np.zeros(self.nwx + 1)
+
+ self.hy_left_0 = np.zeros(self.nwy + 1)
+ self.hy_left_1 = np.zeros(self.nwy + 1)
+
+ self.hy_right_0 = np.zeros(self.nwy + 1)
+ self.hy_right_1 = np.zeros(self.nwy + 1)
+
+ # E precursors
+ self.ez_front_0 = np.zeros(self.nwx + 1)
+ self.ez_front_1 = np.zeros(self.nwx + 1)
+
+ self.ez_back_0 = np.zeros(self.nwx + 1)
+ self.ez_back_1 = np.zeros(self.nwx + 1)
+
+ self.ez_left_0 = np.zeros(self.nwy + 1)
+ self.ez_left_1 = np.zeros(self.nwy + 1)
+
+ self.ez_right_0 = np.zeros(self.nwy + 1)
+ self.ez_right_1 = np.zeros(self.nwy + 1)
+
+ def _initialize_field_names(self):
+ self.fn_m = ['hy_left', 'hy_right', 'hx_back', 'hx_front']
+ self.fn_e = ['ez_left', 'ez_right', 'ez_back', 'ez_front']
+
+ def interpolate_across_sub_cells(self, y, n0, n1):
+ n = n1 - n0
+ x = np.arange(0, y.shape[0], 1.0)
+ x_new = np.linspace(1, y.shape[0] - 2, n * self.ratio + 1)
+ f = interpolate.interp1d(x, y, kind='linear')
+ y_new = f(x_new)
+ """plt.plot(x, y, 'x', x_new, y_new, '.')
+ plt.show()
+ sys.exit()"""
+ return y_new
+
+ def update_electric(self):
+ # f = np.sin(np.arange(0, 13))
+
+ # line of Ez nodes along the left edge of the IS
+ self.ez_left_0 = np.copy(self.ez_left_1)
+ # interpolate nodes between two Ez nodes 1 cell beyond/infront the IS
+ f = self.Ez[self.i0, self.j0 - 1:self.j1 + 2, 0]
+ self.ez_left_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+ self.ez_right_0 = np.copy(self.ez_right_1)
+ f = self.Ez[self.i1, self.j0 - 1:self.j1 + 2, 0]
+ self.ez_right_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+ self.ez_front_0 = np.copy(self.ez_front_1)
+ f = self.Ez[self.i0 - 1:self.i1 + 2, self.j0, 0]
+ self.ez_front_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+ self.ez_back_0 = np.copy(self.ez_back_1)
+ f = self.Ez[self.i0 - 1:self.i1 + 2, self.j1, 0]
+ self.ez_back_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+ def update_hy_left_1(self):
+
+ self.hy_left_0 = np.copy(self.hy_left_1)
+ c1, c2 = calculate_weighting_coefficients(self.l_weight, self.ratio)
+ l = self.Hy[self.i0 - 1, self.j0 - 1:self.j1 + 2, 0]
+ r = self.Hy[self.i0, self.j0 - 1:self.j1 + 2, 0]
+
+ # Tranverse interpolated hz
+ hz_t_i = c1 * l + c2 * r
+ # Tangential interpolated hz
+ self.hy_left_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+ def update_hy_right_1(self):
+
+ self.hy_right_0 = np.copy(self.hy_right_1)
+ c1, c2 = calculate_weighting_coefficients(self.r_weight, self.ratio)
+ l = self.Hy[self.i1 - 1, self.j0 - 1:self.j1 + 2, 0]
+ r = self.Hy[self.i1, self.j0 - 1:self.j1 + 2, 0]
+ hz_t_i = c1 * l + c2 * r
+ self.hy_right_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+ def update_hx_back_1(self):
+ self.hx_back_0 = np.copy(self.hx_back_1)
+ c1, c2 = calculate_weighting_coefficients(self.r_weight, self.ratio)
+ l = self.Hx[self.i0 - 1:self.i1 + 2, self.j1 - 1, 0]
+ r = self.Hx[self.i0 - 1:self.i1 + 2, self.j1, 0]
+ hz_t_i = c1 * l + c2 * r
+ self.hx_back_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+ def update_hx_front_1(self):
+
+ self.hx_front_0 = np.copy(self.hx_front_1)
+ c1, c2 = calculate_weighting_coefficients(self.l_weight, self.ratio)
+ l = self.Hx[self.i0 - 1:self.i1 + 2, self.j0 - 1, 0]
+ r = self.Hx[self.i0 - 1:self.i1 + 2, self.j0, 0]
+ hz_t_i = c1 * l + c2 * r
+ self.hx_front_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+ def update_magnetic(self):
+ self.update_hy_left_1()
+ self.update_hy_right_1()
+ self.update_hx_back_1()
+ self.update_hx_front_1()
+
+
+class PrecursorNodes2d(PrecusorNodes2dBase):
+
+ def __init__(self, fdtd_grid, sub_grid):
+ super().__init__(fdtd_grid, sub_grid)
+
+ def _initialize_fields(self):
+ # E precursors
+ self.ex_front_0 = np.zeros(self.nwx)
+ self.ex_front_1 = np.zeros(self.nwx)
+
+ self.ex_back_0 = np.zeros(self.nwx)
+ self.ex_back_1 = np.zeros(self.nwx)
+
+ self.ey_left_0 = np.zeros(self.nwy)
+ self.ey_left_1 = np.zeros(self.nwy)
+
+ self.ey_right_0 = np.zeros(self.nwy)
+ self.ey_right_1 = np.zeros(self.nwy)
+
+ # H precursors
+ self.hz_front_0 = np.zeros(self.nwx)
+ self.hz_front_1 = np.zeros(self.nwx)
+
+ self.hz_back_0 = np.zeros(self.nwx)
+ self.hz_back_1 = np.zeros(self.nwx)
+
+ self.hz_left_0 = np.zeros(self.nwy)
+ self.hz_left_1 = np.zeros(self.nwy)
+
+ self.hz_right_0 = np.zeros(self.nwy)
+ self.hz_right_1 = np.zeros(self.nwy)
+
+ def _initialize_field_names(self):
+ # Field names
+ self.fn_m = ['hz_left', 'hz_right', 'hz_back', 'hz_front']
+ self.fn_e = ['ey_left', 'ey_right', 'ex_back', 'ex_front']
+
+ def update_hz_left_1(self):
+ self.hz_left_0 = np.copy(self.hz_left_1)
+ c1, c2 = calculate_weighting_coefficients(1, self.ratio)
+ l = self.Hz[self.i0 - 1, self.j0 - 1:self.j1 + 1, 1]
+ r = self.Hz[self.i0, self.j0 - 1:self.j1 + 1, 1]
+
+ # Tranverse interpolated hz
+ hz_t_i = c1 * l + c2 * r
+ # Tangential interpolated hz
+ self.hz_left_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+ def update_hz_right_1(self):
+
+ self.hz_right_0 = np.copy(self.hz_right_1)
+ c1, c2 = calculate_weighting_coefficients(2, self.ratio)
+ l = self.Hz[self.i1 - 1, self.j0 - 1:self.j1 + 1, 1]
+ r = self.Hz[self.i1, self.j0 - 1:self.j1 + 1, 1]
+ hz_t_i = c1 * l + c2 * r
+ self.hz_right_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+ def update_hz_back_1(self):
+ self.hz_back_0 = np.copy(self.hz_back_1)
+ c1, c2 = calculate_weighting_coefficients(2, self.ratio)
+ l = self.Hz[self.i0 - 1:self.i1 + 1, self.j1 - 1, 1]
+ r = self.Hz[self.i0 - 1:self.i1 + 1, self.j1, 1]
+ hz_t_i = c1 * l + c2 * r
+ self.hz_back_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+ def update_hz_front_1(self):
+ self.hz_front_0 = np.copy(self.hz_front_1)
+ c1, c2 = calculate_weighting_coefficients(1, self.ratio)
+ l = self.Hz[self.i0 - 1:self.i1 + 1, self.j0 - 1, 1]
+ r = self.Hz[self.i0 - 1:self.i1 + 1, self.j0, 1]
+ hz_t_i = c1 * l + c2 * r
+ self.hz_front_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+ def update_magnetic(self):
+ self.update_hz_left_1()
+ self.update_hz_right_1()
+ self.update_hz_back_1()
+ self.update_hz_front_1()
+
+ def update_electric(self):
+
+ """Function to calculate the precursor nodes at the next main
+ grid time-step
+ """
+
+ # LEFT BOUNDARY
+ # Save the last precursor node values - these will be used interpolate
+ # field values at sub grid time step values.
+ self.ey_left_0 = np.copy(self.ey_left_1)
+ # line of Ex nodes along the left edge of the IS
+ f = self.Ey[self.i0, self.j0 - 1:self.j1 + 1, 1]
+ self.ey_left_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+ # RIGHT BOUNDARY
+ self.ey_right_0 = np.copy(self.ey_right_1)
+ f = self.Ey[self.i1, self.j0 - 1:self.j1 + 1, 1]
+ self.ey_right_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+ # BACK BOUNDARY
+ self.ex_back_0 = np.copy(self.ex_back_1)
+ f = self.Ex[self.i0 - 1:self.i1 + 1, self.j1, 1]
+ self.ex_back_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+ # FRONT BOUNDARY
+ self.ex_front_0 = np.copy(self.ex_front_1)
+ f = self.Ex[self.i0 - 1:self.i1 + 1, self.j0, 1]
+ self.ex_front_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+ def interpolate_across_sub_cells(self, y, n0, n1):
+ n = n1 - n0
+ x = np.arange(0.5, y.shape[0], 1.0)
+ x_new = np.linspace(1 + self.d, (1 + n) - self.d, n * self.ratio)
+ f = interpolate.interp1d(x, y, kind='linear')
+ y_new = f(x_new)
+
+ return y_new
+
+
+class PrecursorNodes(PrecusorNodes2dBase):
+
+ def __init__(self, fdtd_grid, sub_grid):
+
+ self.nwz = sub_grid.nwz
+ super().__init__(fdtd_grid, sub_grid)
+ self.initialize_magnetic_slices_array()
+ self.initialize_electric_slices_array()
+
+ def build_coefficient_matrix(self, update_coefficients, lookup_id, inc_field, name):
+ """Function builds a 2d matrix of update coefficients for each face and field.
+ This allows the is nodes to be updated via slicing rather than iterating which is much faster
+ """
+ pass
def _initialize_fields(self):
@@ -157,124 +433,12 @@ class PrecusorNodesBase:
'ex_bottom', 'ey_bottom'
]
- def interpolate_magnetic_in_time(self, m):
- self.weight_pre_and_current_fields(m, self.fn_m)
-
- def interpolate_electric_in_time(self, m):
- self.weight_pre_and_current_fields(m, self.fn_e)
-
- def weight_pre_and_current_fields(self, m, field_names):
- c1, c2 = calculate_weighting_coefficients(m, self.ratio)
-
- for f in field_names:
- try:
- val = c1 * getattr(self, f + '_0') + c2 * getattr(self, f + '_1')
- except ValueError:
- print(self.ex_front_0.shape)
- print(self.ex_front_1.shape)
- raise Exception(f)
- setattr(self, f, val)
-
- def calc_exact_field(self, field_names):
- """
- Function to set the fields used in update calculations to the
- values at the current main time step.
- i.e. ey_left = copy.ey_left_1
- """
- for f in field_names:
- val = np.copy(getattr(self, f + '_1'))
- setattr(self, f, val)
-
def update_previous_timestep_fields(self, field_names):
for fn in field_names:
val = getattr(self, fn + '_1')
val_c = np.copy(val)
setattr(self, fn + '_0', val_c)
- def calc_exact_magnetic_in_time(self):
- self.calc_exact_field(self.fn_m)
-
- def calc_exact_electric_in_time(self):
- self.calc_exact_field(self.fn_e)
-
- def create_interpolated_coords(self, mid, field):
-
- n_x = field.shape[0]
- n_y = field.shape[1]
-
- if mid:
- x = np.arange(0.5, n_x, 1.0)
- z = np.arange(0, n_y, 1.0)
-
- # Coordinates that require interpolated values
- x_sg = np.linspace(self.d, n_x - self.d, n_x * self.ratio)
- z_sg = np.linspace(0, n_y - 1, (n_y - 1) * self.ratio + 1)
-
- else:
- x = np.arange(0, n_x, 1.0)
- z = np.arange(0.5, n_y, 1.0)
-
- # Coordinates that require interpolated values
- x_sg = np.linspace(0, n_x - 1, (n_x - 1) * self.ratio + 1)
- z_sg = np.linspace(self.d, n_y - self.d, n_y * self.ratio)
-
- return (x, z, x_sg, z_sg)
-
- def interpolate_to_sub_grid(self, field, coords):
- x, z, x_sg, z_sg = coords
- interp_f = interpolate.RectBivariateSpline(x,
- z,
- field,
- kx=self.interpolation,
- ky=self.interpolation)
- f_i = interp_f(x_sg, z_sg)
- return f_i
-
- def update_magnetic(self):
-
- # Copy previous time step magnetic field values to the previous
- # time step variables
- self.update_previous_timestep_fields(self.fn_m)
-
- for obj in self.magnetic_slices:
-
- # Grab the main grid fields used to interpolate across the IS
- # f = self.Hi[slice]
- f_1, f_2 = self.f1_f2_magnetic(obj)
- w = obj[-2]
- c1, c2 = calculate_weighting_coefficients(w, self.ratio)
- # transverse interpolated h field
- f_t = c1 * f_1 + c2 * f_2
-
- # interpolate over a fine grid
- f_i = self.interpolate_to_sub_grid(f_t, obj[1])
-
- if f_i == f_t:
- raise ValueError
- # discard the outer nodes only required for interpolation
- f = f_i[self.ratio:-self.ratio, self.ratio:-self.ratio]
- setattr(self, obj[0], f)
-
- def update_electric(self):
-
- self.update_previous_timestep_fields(self.fn_e)
-
- for obj in self.electric_slices:
- 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]
- setattr(self, obj[0], f)
-
-
-class PrecursorNodesFiltered(PrecusorNodesBase):
-
- def __init__(self, fdtd_grid, sub_grid):
-
- super().__init__(fdtd_grid, sub_grid)
- 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
@@ -283,99 +447,59 @@ class PrecursorNodesFiltered(PrecusorNodesBase):
# Extend the surface so that the outer fields can be interpolated
# more accurately
- i0 = self.i0 - 1
- j0 = self.j0 - 1
- k0 = self.k0 - 1
- i1 = self.i1 + 1
- j1 = self.j1 + 1
- k1 = self.k1 + 1
+ #i0 = self.i0 - 1
+ #j0 = self.j0 - 1
+ #k0 = self.k0 - 1
+ #i1 = self.i1 + 1
+ #j1 = self.j1 + 1
+ #k1 = self.k1 + 1
+
+ # not extended
+ i0 = self.i0
+ j0 = self.j0
+ k0 = self.k0
+ i1 = self.i1
+ j1 = self.j1
+ k1 = self.k1
slices = [
['hy_left_1', False,
- (self.i0 - 2, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
(self.i0 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- (self.i0 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- self.l_weight,
- self.Hy],
+ (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)), self.Hy],
['hy_right_1', False,
- (self.i1 - 2, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
(self.i1 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- (self.i1 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- self.r_weight,
- self.Hy],
+ (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)), self.Hy],
['hz_left_1', True,
- (self.i0 - 2, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
(self.i0 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- (self.i0 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- self.l_weight,
- self.Hz],
+ (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)), self.Hz],
['hz_right_1', True,
- (self.i1 - 2, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
(self.i1 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- (self.i1 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- self.r_weight,
- self.Hz],
+ (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)), self.Hz],
['hx_front_1', False,
- (slice(i0, i1 + 1, 1), self.j0 - 2, slice(k0, k1, 1)),
(slice(i0, i1 + 1, 1), self.j0 - 1, slice(k0, k1, 1)),
- (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)),
- (slice(i0, i1 + 1, 1), self.j0 + 1, slice(k0, k1, 1)),
- self.l_weight,
- self.Hx],
+ (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)), self.Hx],
['hx_back_1', False,
- (slice(i0, i1 + 1, 1), self.j1 - 2, slice(k0, k1, 1)),
(slice(i0, i1 + 1, 1), self.j1 - 1, slice(k0, k1, 1)),
- (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)),
- (slice(i0, i1 + 1, 1), self.j1 + 1, slice(k0, k1, 1)),
- self.r_weight,
- self.Hx],
+ (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)), self.Hx],
['hz_front_1', True,
- (slice(i0, i1, 1), self.j0 - 2, slice(k0, k1 + 1, 1)),
(slice(i0, i1, 1), self.j0 - 1, slice(k0, k1 + 1, 1)),
- (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)),
- (slice(i0, i1, 1), self.j0 + 1, slice(k0, k1 + 1, 1)),
- self.l_weight,
- self.Hz],
+ (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)), self.Hz],
['hz_back_1', True,
- (slice(i0, i1, 1), self.j1 - 2, slice(k0, k1 + 1, 1)),
(slice(i0, i1, 1), self.j1 - 1, slice(k0, k1 + 1, 1)),
- (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)),
- (slice(i0, i1, 1), self.j1 + 1, slice(k0, k1 + 1, 1)),
- self.r_weight,
- self.Hz],
+ (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)), self.Hz],
['hx_bottom_1', False,
# check these indexes
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 2),
(slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 1),
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0),
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 + 1),
- self.l_weight,
- self.Hx],
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0), self.Hx],
['hx_top_1', False,
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 2),
(slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 1),
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1),
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 + 1),
- self.r_weight,
- self.Hx],
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1), self.Hx],
['hy_bottom_1', True,
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 2),
(slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 1),
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0),
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 + 1),
- self.l_weight,
- self.Hy],
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0), self.Hy],
['hy_top_1', True,
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 2),
(slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 1),
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1),
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 + 1),
- self.r_weight,
- self.Hy]
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1), self.Hy]
]
for obj in slices:
@@ -387,208 +511,20 @@ class PrecursorNodesFiltered(PrecusorNodesBase):
def initialize_electric_slices_array(self):
# Extend the region sliced from the main grid by 1 cell.
# this allows more accurate interpolation for the outernodes
- i0 = self.i0 - 1
- j0 = self.j0 - 1
- k0 = self.k0 - 1
- i1 = self.i1 + 1
- j1 = self.j1 + 1
- k1 = self.k1 + 1
+ #i0 = self.i0 - 1
+ #j0 = self.j0 - 1
+ #k0 = self.k0 - 1
+ #i1 = self.i1 + 1
+ #j1 = self.j1 + 1
+ #k1 = self.k1 + 1
- # Spatially interpolate nodes
- slices = [
- ['ex_front_1', True,
- (slice(i0, i1, 1), self.j0 - 1, slice(k0, k1 + 1, 1)),
- (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)),
- (slice(i0, i1, 1), self.j0 + 1, slice(k0, k1 + 1, 1)),
- self.Ex],
- ['ex_back_1', True,
- (slice(i0, i1, 1), self.j1 - 1, slice(k0, k1 + 1, 1)),
- (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)),
- (slice(i0, i1, 1), self.j1 + 1, slice(k0, k1 + 1, 1)),
- self.Ex],
- ['ez_front_1', False,
- (slice(i0, i1 + 1, 1), self.j0 - 1, slice(k0, k1, 1)),
- (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)),
- (slice(i0, i1 + 1, 1), self.j0 + 1, slice(k0, k1, 1)),
- self.Ez],
- ['ez_back_1', False,
- (slice(i0, i1 + 1, 1), self.j1 - 1, slice(k0, k1, 1)),
- (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)),
- (slice(i0, i1 + 1, 1), self.j1 + 1, slice(k0, k1, 1)),
- self.Ez],
- ['ey_left_1', True,
- (self.i0 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- (self.i0 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- self.Ey],
- ['ey_right_1', True,
- (self.i1 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- (self.i1 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- self.Ey],
- ['ez_left_1', False,
- (self.i0 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- (self.i0 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- self.Ez],
- ['ez_right_1', False,
- (self.i1 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- (self.i1 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- self.Ez],
-
- ['ex_bottom_1', True,
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 1),
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0),
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 + 1),
- self.Ex],
- ['ex_top_1', True,
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 1),
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1),
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 + 1),
- self.Ex],
- ['ey_bottom_1', False,
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 1),
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0),
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 + 1),
- self.Ey],
- ['ey_top_1', False,
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 1),
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1),
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 + 1),
- self.Ey]
- ]
-
- for obj in slices:
- sliced_field = obj[-1][obj[2]]
- obj[1] = self.create_interpolated_coords(obj[1], sliced_field)
-
- self.electric_slices = slices
-
- def f1_f2_magnetic(self, obj):
- f_u_1 = obj[-1][obj[2]]
- f_u_2 = obj[-1][obj[3]]
- f_u_3 = obj[-1][obj[4]]
- f_u_4 = obj[-1][obj[5]]
-
- f_1 = 0.25 * f_u_1 + 0.5 * f_u_2 + 0.25 * f_u_3
- f_2 = 0.25 * f_u_2 + 0.5 * f_u_3 + 0.25 * f_u_4
-
- return f_1, f_2
-
- def f_m_electric(self, obj):
-
- f_u_1 = obj[-1][obj[2]]
- f_u_2 = obj[-1][obj[3]]
- f_u_3 = obj[-1][obj[4]]
-
- f_m = 0.25 * f_u_1 + 0.5 * f_u_2 + 0.25 * f_u_3
-
- return f_m
-
-
-class PrecursorNodes(PrecusorNodesBase):
-
- def __init__(self, fdtd_grid, sub_grid):
-
- super().__init__(fdtd_grid, sub_grid)
- 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
- # to obtain the 2 2d array of H nodes required for interpolation
- # across the IS boundary for each h field on each face of the subgrid
-
- # Extend the surface so that the outer fields can be interpolated
- # more accurately
- i0 = self.i0 - 1
- j0 = self.j0 - 1
- k0 = self.k0 - 1
- i1 = self.i1 + 1
- j1 = self.j1 + 1
- k1 = self.k1 + 1
-
- slices = [
- ['hy_left_1', False,
- (self.i0 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- self.l_weight,
- self.Hy],
- ['hy_right_1', False,
- (self.i1 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
- self.r_weight,
- self.Hy],
- ['hz_left_1', True,
- (self.i0 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- self.l_weight,
- self.Hz],
- ['hz_right_1', True,
- (self.i1 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
- self.r_weight,
- self.Hz],
- ['hx_front_1', False,
- (slice(i0, i1 + 1, 1), self.j0 - 1, slice(k0, k1, 1)),
- (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)),
- self.l_weight,
- self.Hx],
- ['hx_back_1', False,
- (slice(i0, i1 + 1, 1), self.j1 - 1, slice(k0, k1, 1)),
- (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)),
- self.r_weight,
- self.Hx],
- ['hz_front_1', True,
- (slice(i0, i1, 1), self.j0 - 1, slice(k0, k1 + 1, 1)),
- (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)),
- self.l_weight,
- self.Hz],
- ['hz_back_1', True,
- (slice(i0, i1, 1), self.j1 - 1, slice(k0, k1 + 1, 1)),
- (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)),
- self.r_weight,
- self.Hz],
- ['hx_bottom_1', False,
- # check these indexes
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 1),
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0),
- self.l_weight,
- self.Hx],
- ['hx_top_1', False,
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 1),
- (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1),
- self.r_weight,
- self.Hx],
- ['hy_bottom_1', True,
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 1),
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0),
- self.l_weight,
- self.Hy],
- ['hy_top_1', True,
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 1),
- (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1),
- self.r_weight,
- self.Hy]
- ]
-
- for obj in slices:
- sliced_field = obj[-1][obj[2]]
- obj[1] = self.create_interpolated_coords(obj[1], sliced_field)
-
- self.magnetic_slices = slices
-
- def initialize_electric_slices_array(self):
- # Extend the region sliced from the main grid by 1 cell.
- # this allows more accurate interpolation for the outernodes
- i0 = self.i0 - 1
- j0 = self.j0 - 1
- k0 = self.k0 - 1
- i1 = self.i1 + 1
- j1 = self.j1 + 1
- k1 = self.k1 + 1
+ # not extended
+ i0 = self.i0
+ j0 = self.j0
+ k0 = self.k0
+ i1 = self.i1
+ j1 = self.j1
+ k1 = self.k1
# Spatially interpolate nodes
slices = [
@@ -614,10 +550,84 @@ class PrecursorNodes(PrecusorNodesBase):
self.electric_slices = slices
- def f1_f2_magnetic(self, obj):
- f_1 = obj[-1][obj[2]]
- f_2 = obj[-1][obj[3]]
- return f_1, f_2
+ def create_interpolated_coords(self, mid, field):
- def f_m_electric(self, obj):
- return obj[-1][obj[2]]
+ n_x = field.shape[0]
+ n_y = field.shape[1]
+
+ if mid:
+ x = np.arange(0.5, n_x, 1.0)
+ z = np.arange(0, n_y, 1.0)
+
+ # Coordinates that require interpolated values
+ x_sg = np.linspace(self.d, n_x - self.d, n_x * self.ratio)
+ z_sg = np.linspace(0, n_y - 1, (n_y - 1) * self.ratio + 1)
+
+ else:
+ x = np.arange(0, n_x, 1.0)
+ z = np.arange(0.5, n_y, 1.0)
+
+ # Coordinates that require interpolated values
+ x_sg = np.linspace(0, n_x - 1, (n_x - 1) * self.ratio + 1)
+ z_sg = np.linspace(self.d, n_y - self.d, n_y * self.ratio)
+
+ return (x, z, x_sg, z_sg)
+
+ def update_magnetic(self):
+
+ # Copy previous time step magnetic field values to the previous
+ # time step variables
+ self.update_previous_timestep_fields(self.fn_m)
+
+ for obj in self.magnetic_slices:
+
+ # Grab the main grid fields used to interpolate across the IS
+ # f = self.Hi[slice]
+ f_1 = obj[-1][obj[2]]
+ f_2 = obj[-1][obj[3]]
+ if ('left' in obj[0] or
+ 'bottom' in obj[0] or
+ 'front' in obj[0]):
+ w = self.l_weight
+ else:
+ w = self.r_weight
+ c1, c2 = calculate_weighting_coefficients(w, self.ratio)
+ # transverse interpolated h field
+ f_t = c1 * f_1 + c2 * f_2
+
+ # interpolate over a fine grid
+ f_i = self.interpolate_to_sub_grid(f_t, obj[1])
+
+ if f_i == f_t:
+ raise ValueError
+
+ # discard the outer nodes only required for interpolation
+ #f = f_i[self.ratio:-self.ratio, self.ratio:-self.ratio]
+ f = f_i
+ setattr(self, obj[0], f)
+
+ def update_electric(self):
+
+ self.update_previous_timestep_fields(self.fn_e)
+
+ for obj in self.electric_slices:
+ f_m = obj[-1][obj[2]]
+ f_i = self.interpolate_to_sub_grid(f_m, obj[1])
+ f = f_i
+ #f = f_i[self.ratio:-self.ratio, self.ratio:-self.ratio]
+ setattr(self, obj[0], f)
+
+ def interpolate_to_sub_grid(self, field, coords):
+ x, z, x_sg, z_sg = coords
+ ex_t = np.transpose(field)
+ f = interpolate.interp2d(x, z, ex_t, kind=self.interpolation)
+ #f = interpolate.RectBivariateSpline(x, z, field)
+ ex_inter_t = f(x_sg, z_sg)
+ ex_inter = np.transpose(ex_inter_t)
+ #ex_inter = ex_inter_t
+
+ return ex_inter
+
+
+class PlaneError(Exception):
+ pass
diff --git a/gprMax/subgrids/precursor_nodes_filtered.py b/gprMax/subgrids/precursor_nodes_filtered.py
new file mode 100644
index 00000000..a1046c76
--- /dev/null
+++ b/gprMax/subgrids/precursor_nodes_filtered.py
@@ -0,0 +1,738 @@
+import numpy as np
+from scipy import interpolate
+import sys
+
+
+def calculate_weighting_coefficients(x1, x):
+ c1 = (x - x1) / x
+ c2 = x1 / x
+ return (c1, c2)
+
+
+class PrecusorNodes2dBase(object):
+
+ def __init__(self, fdtd_grid, sub_grid):
+ self.G = fdtd_grid
+ self.ratio = sub_grid.ratio
+ self.nwx = sub_grid.nwx
+ self.nwy = sub_grid.nwy
+ self.sub_grid = sub_grid
+ self.interpolation = sub_grid.interpolation
+
+ self.Hx = fdtd_grid.Hx
+ self.Hy = fdtd_grid.Hy
+ self.Hz = fdtd_grid.Hz
+ self.Ex = fdtd_grid.Ex
+ self.Ey = fdtd_grid.Ey
+ self.Ez = fdtd_grid.Ez
+
+ # Main grid indices of subgrids
+ self.i0 = sub_grid.i0
+ self.j0 = sub_grid.j0
+ self.k0 = sub_grid.k0
+ self.i1 = sub_grid.i1
+ self.j1 = sub_grid.j1
+ self.k1 = sub_grid.k1
+
+ # dl / 2 sub cell
+ self.d = 1 / (2 * self.ratio)
+
+ self._initialize_fields()
+ self._initialize_field_names()
+
+ self.l_weight = self.ratio // 2
+ self.r_weight = self.ratio - self.l_weight
+ #self.l_weight = 1
+ #self.r_weight = 2
+
+ def _initialize_fields(self):
+ print('dont call me')
+ sys.exit()
+ pass
+
+ def _initialize_field_names(self):
+ print('dont call me')
+ sys.exit()
+ pass
+
+ def interpolate_magnetic_in_time(self, m):
+ self.weight_pre_and_current_fields(m, self.fn_m)
+
+ def interpolate_electric_in_time(self, m):
+ self.weight_pre_and_current_fields(m, self.fn_e)
+
+ def weight_pre_and_current_fields(self, m, field_names):
+ c1, c2 = calculate_weighting_coefficients(m, self.ratio)
+
+ for f in field_names:
+ try:
+ val = c1 * getattr(self, f + '_0') + c2 * getattr(self, f + '_1')
+ except ValueError:
+ print(self.ex_front_0.shape)
+ print(self.ex_front_1.shape)
+ raise Exception(f)
+ setattr(self, f, val)
+
+ def calc_exact_field(self, field_names):
+ """
+ Function to set the fields used in update calculations to the
+ values at the current main time step.
+ i.e. ey_left = copy.ey_left_1
+ """
+ for f in field_names:
+ val = np.copy(getattr(self, f + '_1'))
+ setattr(self, f, val)
+
+ def calc_exact_magnetic_in_time(self):
+ self.calc_exact_field(self.fn_m)
+
+ def calc_exact_electric_in_time(self):
+ self.calc_exact_field(self.fn_e)
+
+
+class PrecursorNodes2dTM(PrecusorNodes2dBase):
+
+ def __init__(self, fdtd_grid, sub_grid):
+ super().__init__(fdtd_grid, sub_grid)
+
+ def _initialize_fields(self):
+
+ # H precursors
+ self.hx_front_0 = np.zeros(self.nwx + 1)
+ self.hx_front_1 = np.zeros(self.nwx + 1)
+
+ self.hx_back_0 = np.zeros(self.nwx + 1)
+ self.hx_back_1 = np.zeros(self.nwx + 1)
+
+ self.hy_left_0 = np.zeros(self.nwy + 1)
+ self.hy_left_1 = np.zeros(self.nwy + 1)
+
+ self.hy_right_0 = np.zeros(self.nwy + 1)
+ self.hy_right_1 = np.zeros(self.nwy + 1)
+
+ # E precursors
+ self.ez_front_0 = np.zeros(self.nwx + 1)
+ self.ez_front_1 = np.zeros(self.nwx + 1)
+
+ self.ez_back_0 = np.zeros(self.nwx + 1)
+ self.ez_back_1 = np.zeros(self.nwx + 1)
+
+ self.ez_left_0 = np.zeros(self.nwy + 1)
+ self.ez_left_1 = np.zeros(self.nwy + 1)
+
+ self.ez_right_0 = np.zeros(self.nwy + 1)
+ self.ez_right_1 = np.zeros(self.nwy + 1)
+
+ def _initialize_field_names(self):
+ self.fn_m = ['hy_left', 'hy_right', 'hx_back', 'hx_front']
+ self.fn_e = ['ez_left', 'ez_right', 'ez_back', 'ez_front']
+
+ def interpolate_across_sub_cells(self, y, n0, n1):
+ n = n1 - n0
+ x = np.arange(0, y.shape[0], 1.0)
+ x_new = np.linspace(1, y.shape[0] - 2, n * self.ratio + 1)
+ f = interpolate.interp1d(x, y, kind='linear')
+ y_new = f(x_new)
+ """plt.plot(x, y, 'x', x_new, y_new, '.')
+ plt.show()
+ sys.exit()"""
+ return y_new
+
+ def update_electric(self):
+ # f = np.sin(np.arange(0, 13))
+
+ # line of Ez nodes along the left edge of the IS
+ self.ez_left_0 = np.copy(self.ez_left_1)
+ # interpolate nodes between two Ez nodes 1 cell beyond/infront the IS
+ f = self.Ez[self.i0, self.j0 - 1:self.j1 + 2, 0]
+ self.ez_left_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+ self.ez_right_0 = np.copy(self.ez_right_1)
+ f = self.Ez[self.i1, self.j0 - 1:self.j1 + 2, 0]
+ self.ez_right_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+ self.ez_front_0 = np.copy(self.ez_front_1)
+ f = self.Ez[self.i0 - 1:self.i1 + 2, self.j0, 0]
+ self.ez_front_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+ self.ez_back_0 = np.copy(self.ez_back_1)
+ f = self.Ez[self.i0 - 1:self.i1 + 2, self.j1, 0]
+ self.ez_back_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+ def update_hy_left_1(self):
+
+ self.hy_left_0 = np.copy(self.hy_left_1)
+ c1, c2 = calculate_weighting_coefficients(self.l_weight, self.ratio)
+ l = self.Hy[self.i0 - 1, self.j0 - 1:self.j1 + 2, 0]
+ r = self.Hy[self.i0, self.j0 - 1:self.j1 + 2, 0]
+
+ # Tranverse interpolated hz
+ hz_t_i = c1 * l + c2 * r
+ # Tangential interpolated hz
+ self.hy_left_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+ def update_hy_right_1(self):
+
+ self.hy_right_0 = np.copy(self.hy_right_1)
+ c1, c2 = calculate_weighting_coefficients(self.r_weight, self.ratio)
+ l = self.Hy[self.i1 - 1, self.j0 - 1:self.j1 + 2, 0]
+ r = self.Hy[self.i1, self.j0 - 1:self.j1 + 2, 0]
+ hz_t_i = c1 * l + c2 * r
+ self.hy_right_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+ def update_hx_back_1(self):
+ self.hx_back_0 = np.copy(self.hx_back_1)
+ c1, c2 = calculate_weighting_coefficients(self.r_weight, self.ratio)
+ l = self.Hx[self.i0 - 1:self.i1 + 2, self.j1 - 1, 0]
+ r = self.Hx[self.i0 - 1:self.i1 + 2, self.j1, 0]
+ hz_t_i = c1 * l + c2 * r
+ self.hx_back_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+ def update_hx_front_1(self):
+
+ self.hx_front_0 = np.copy(self.hx_front_1)
+ c1, c2 = calculate_weighting_coefficients(self.l_weight, self.ratio)
+ l = self.Hx[self.i0 - 1:self.i1 + 2, self.j0 - 1, 0]
+ r = self.Hx[self.i0 - 1:self.i1 + 2, self.j0, 0]
+ hz_t_i = c1 * l + c2 * r
+ self.hx_front_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+ def update_magnetic(self):
+ self.update_hy_left_1()
+ self.update_hy_right_1()
+ self.update_hx_back_1()
+ self.update_hx_front_1()
+
+
+class PrecursorNodes2d(PrecusorNodes2dBase):
+
+ def __init__(self, fdtd_grid, sub_grid):
+ super().__init__(fdtd_grid, sub_grid)
+
+ def _initialize_fields(self):
+ # E precursors
+ self.ex_front_0 = np.zeros(self.nwx)
+ self.ex_front_1 = np.zeros(self.nwx)
+
+ self.ex_back_0 = np.zeros(self.nwx)
+ self.ex_back_1 = np.zeros(self.nwx)
+
+ self.ey_left_0 = np.zeros(self.nwy)
+ self.ey_left_1 = np.zeros(self.nwy)
+
+ self.ey_right_0 = np.zeros(self.nwy)
+ self.ey_right_1 = np.zeros(self.nwy)
+
+ # H precursors
+ self.hz_front_0 = np.zeros(self.nwx)
+ self.hz_front_1 = np.zeros(self.nwx)
+
+ self.hz_back_0 = np.zeros(self.nwx)
+ self.hz_back_1 = np.zeros(self.nwx)
+
+ self.hz_left_0 = np.zeros(self.nwy)
+ self.hz_left_1 = np.zeros(self.nwy)
+
+ self.hz_right_0 = np.zeros(self.nwy)
+ self.hz_right_1 = np.zeros(self.nwy)
+
+ def _initialize_field_names(self):
+ # Field names
+ self.fn_m = ['hz_left', 'hz_right', 'hz_back', 'hz_front']
+ self.fn_e = ['ey_left', 'ey_right', 'ex_back', 'ex_front']
+
+ def update_hz_left_1(self):
+ self.hz_left_0 = np.copy(self.hz_left_1)
+ c1, c2 = calculate_weighting_coefficients(1, self.ratio)
+ l = self.Hz[self.i0 - 1, self.j0 - 1:self.j1 + 1, 1]
+ r = self.Hz[self.i0, self.j0 - 1:self.j1 + 1, 1]
+
+ # Tranverse interpolated hz
+ hz_t_i = c1 * l + c2 * r
+ # Tangential interpolated hz
+ self.hz_left_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+ def update_hz_right_1(self):
+
+ self.hz_right_0 = np.copy(self.hz_right_1)
+ c1, c2 = calculate_weighting_coefficients(2, self.ratio)
+ l = self.Hz[self.i1 - 1, self.j0 - 1:self.j1 + 1, 1]
+ r = self.Hz[self.i1, self.j0 - 1:self.j1 + 1, 1]
+ hz_t_i = c1 * l + c2 * r
+ self.hz_right_1 = self.interpolate_across_sub_cells(hz_t_i, self.j0, self.j1)
+
+ def update_hz_back_1(self):
+ self.hz_back_0 = np.copy(self.hz_back_1)
+ c1, c2 = calculate_weighting_coefficients(2, self.ratio)
+ l = self.Hz[self.i0 - 1:self.i1 + 1, self.j1 - 1, 1]
+ r = self.Hz[self.i0 - 1:self.i1 + 1, self.j1, 1]
+ hz_t_i = c1 * l + c2 * r
+ self.hz_back_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+ def update_hz_front_1(self):
+ self.hz_front_0 = np.copy(self.hz_front_1)
+ c1, c2 = calculate_weighting_coefficients(1, self.ratio)
+ l = self.Hz[self.i0 - 1:self.i1 + 1, self.j0 - 1, 1]
+ r = self.Hz[self.i0 - 1:self.i1 + 1, self.j0, 1]
+ hz_t_i = c1 * l + c2 * r
+ self.hz_front_1 = self.interpolate_across_sub_cells(hz_t_i, self.i0, self.i1)
+
+ def update_magnetic(self):
+ self.update_hz_left_1()
+ self.update_hz_right_1()
+ self.update_hz_back_1()
+ self.update_hz_front_1()
+
+ def update_electric(self):
+
+ """Function to calculate the precursor nodes at the next main
+ grid time-step
+ """
+
+ # LEFT BOUNDARY
+ # Save the last precursor node values - these will be used interpolate
+ # field values at sub grid time step values.
+ self.ey_left_0 = np.copy(self.ey_left_1)
+ # line of Ex nodes along the left edge of the IS
+ f = self.Ey[self.i0, self.j0 - 1:self.j1 + 1, 1]
+ self.ey_left_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+ # RIGHT BOUNDARY
+ self.ey_right_0 = np.copy(self.ey_right_1)
+ f = self.Ey[self.i1, self.j0 - 1:self.j1 + 1, 1]
+ self.ey_right_1 = self.interpolate_across_sub_cells(f, self.j0, self.j1)
+
+ # BACK BOUNDARY
+ self.ex_back_0 = np.copy(self.ex_back_1)
+ f = self.Ex[self.i0 - 1:self.i1 + 1, self.j1, 1]
+ self.ex_back_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+ # FRONT BOUNDARY
+ self.ex_front_0 = np.copy(self.ex_front_1)
+ f = self.Ex[self.i0 - 1:self.i1 + 1, self.j0, 1]
+ self.ex_front_1 = self.interpolate_across_sub_cells(f, self.i0, self.i1)
+
+ def interpolate_across_sub_cells(self, y, n0, n1):
+ n = n1 - n0
+ x = np.arange(0.5, y.shape[0], 1.0)
+ x_new = np.linspace(1 + self.d, (1 + n) - self.d, n * self.ratio)
+ f = interpolate.interp1d(x, y, kind='linear')
+ y_new = f(x_new)
+
+ return y_new
+
+
+class PrecursorNodes(PrecusorNodes2dBase):
+
+ def __init__(self, fdtd_grid, sub_grid):
+
+ self.nwz = sub_grid.nwz
+ super().__init__(fdtd_grid, sub_grid)
+ self.initialize_magnetic_slices_array()
+ self.initialize_electric_slices_array()
+
+ def build_coefficient_matrix(self, update_coefficients, lookup_id, inc_field, name):
+ """Function builds a 2d matrix of update coefficients for each face and field.
+ This allows the is nodes to be updated via slicing rather than iterating which is much faster
+ """
+ pass
+
+ def _initialize_fields(self):
+
+ # Initialise the precursor arrays
+
+ # The precursors are divided up into the 6. Each represent 1
+ # face of a huygens cube surface. We represent each face as a 2d array
+ # containing a field in a particular direction.
+
+ # _1 are the fields at the current main grid timestep
+ # _0 are the fields at the previous main grid timestep
+ # We store both fields so we can do different interpolations between
+ # them on the fly.
+
+ # Front face
+ self.ex_front_1 = np.zeros((self.nwx, self.nwz + 1))
+ self.ex_front_0 = np.copy(self.ex_front_1)
+ self.ez_front_1 = np.zeros((self.nwx + 1, self.nwz))
+ self.ez_front_0 = np.copy(self.ez_front_1)
+
+ # The same as the opposite face
+ self.ex_back_1 = np.copy(self.ex_front_1)
+ self.ex_back_0 = np.copy(self.ex_front_1)
+ self.ez_back_1 = np.copy(self.ez_front_1)
+ self.ez_back_0 = np.copy(self.ez_front_1)
+
+ self.ey_left_1 = np.zeros((self.nwy, self.nwz + 1))
+ self.ey_left_0 = np.copy(self.ey_left_1)
+ self.ez_left_1 = np.zeros((self.nwy + 1, self.nwz))
+ self.ez_left_0 = np.copy(self.ez_left_1)
+
+ self.ey_right_1 = np.copy(self.ey_left_1)
+ self.ey_right_0 = np.copy(self.ey_left_1)
+ self.ez_right_1 = np.copy(self.ez_left_1)
+ self.ez_right_0 = np.copy(self.ez_left_1)
+
+ self.ex_bottom_1 = np.zeros((self.nwx, self.nwy + 1))
+ self.ex_bottom_0 = np.copy(self.ex_bottom_1)
+ self.ey_bottom_1 = np.zeros((self.nwx + 1, self.nwy))
+ self.ey_bottom_0 = np.copy(self.ey_bottom_1)
+
+ self.ex_top_1 = np.copy(self.ex_bottom_1)
+ self.ex_top_0 = np.copy(self.ex_bottom_1)
+ self.ey_top_1 = np.copy(self.ey_bottom_1)
+ self.ey_top_0 = np.copy(self.ey_bottom_1)
+
+ # Initialize the H precursor fields
+ self.hx_front_1 = np.copy(self.ez_front_1)
+ self.hx_front_0 = np.copy(self.ez_front_1)
+ self.hz_front_1 = np.copy(self.ex_front_1)
+ self.hz_front_0 = np.copy(self.ex_front_1)
+
+ self.hx_back_1 = np.copy(self.hx_front_1)
+ self.hx_back_0 = np.copy(self.hx_front_1)
+ self.hz_back_1 = np.copy(self.hz_front_1)
+ self.hz_back_0 = np.copy(self.hz_front_1)
+
+ self.hy_left_1 = np.copy(self.ez_left_1)
+ self.hy_left_0 = np.copy(self.ez_left_1)
+ self.hz_left_1 = np.copy(self.ey_left_1)
+ self.hz_left_0 = np.copy(self.ey_left_1)
+
+ self.hy_right_1 = np.copy(self.hy_left_1)
+ self.hy_right_0 = np.copy(self.hy_left_1)
+ self.hz_right_1 = np.copy(self.hz_left_1)
+ self.hz_right_0 = np.copy(self.hz_left_1)
+
+ self.hx_top_1 = np.copy(self.ey_top_1)
+ self.hx_top_0 = np.copy(self.ey_top_1)
+ self.hy_top_1 = np.copy(self.ex_top_1)
+ self.hy_top_0 = np.copy(self.ex_top_1)
+
+ self.hx_bottom_1 = np.copy(self.hx_top_1)
+ self.hx_bottom_0 = np.copy(self.hx_top_1)
+ self.hy_bottom_1 = np.copy(self.hy_top_1)
+ self.hy_bottom_0 = np.copy(self.hy_top_1)
+
+ def _initialize_field_names(self):
+
+ self.fn_m = [
+ 'hx_front', 'hz_front',
+ 'hx_back', 'hz_back',
+ 'hy_left', 'hz_left',
+ 'hy_right', 'hz_right',
+ 'hx_top', 'hy_top',
+ 'hx_bottom', 'hy_bottom'
+ ]
+
+ self.fn_e = [
+ 'ex_front', 'ez_front',
+ 'ex_back', 'ez_back',
+ 'ey_left', 'ez_left',
+ 'ey_right', 'ez_right',
+ 'ex_top', 'ey_top',
+ 'ex_bottom', 'ey_bottom'
+ ]
+
+ def update_previous_timestep_fields(self, field_names):
+ for fn in field_names:
+ val = getattr(self, fn + '_1')
+ val_c = np.copy(val)
+ setattr(self, fn + '_0', val_c)
+
+ def initialize_magnetic_slices_array(self):
+
+ # Array contains the indices at which the main grid should be sliced
+ # to obtain the 2 2d array of H nodes required for interpolation
+ # across the IS boundary for each h field on each face of the subgrid
+
+ # Extend the surface so that the outer fields can be interpolated
+ # more accurately
+ #i0 = self.i0 - 1
+ #j0 = self.j0 - 1
+ #k0 = self.k0 - 1
+ #i1 = self.i1 + 1
+ #j1 = self.j1 + 1
+ #k1 = self.k1 + 1
+
+ # not extended
+ i0 = self.i0
+ j0 = self.j0
+ k0 = self.k0
+ i1 = self.i1
+ j1 = self.j1
+ k1 = self.k1
+
+ slices = [
+ ['hy_left_1', False,
+ (self.i0 - 2, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+ (self.i0 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+ (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+ (self.i0 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)), self.Hy],
+ ['hy_right_1', False,
+ (self.i1 - 2, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+ (self.i1 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+ (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+ (self.i1 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)), self.Hy],
+ ['hz_left_1', True,
+ (self.i0 - 2, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+ (self.i0 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+ (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+ (self.i0 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)), self.Hz],
+ ['hz_right_1', True,
+ (self.i1 - 2, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+ (self.i1 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+ (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+ (self.i1 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)), self.Hz],
+ ['hx_front_1', False,
+ (slice(i0, i1 + 1, 1), self.j0 - 2, slice(k0, k1, 1)),
+ (slice(i0, i1 + 1, 1), self.j0 - 1, slice(k0, k1, 1)),
+ (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)),
+ (slice(i0, i1 + 1, 1), self.j0 + 1, slice(k0, k1, 1)), self.Hx],
+ ['hx_back_1', False,
+ (slice(i0, i1 + 1, 1), self.j1 - 2, slice(k0, k1, 1)),
+ (slice(i0, i1 + 1, 1), self.j1 - 1, slice(k0, k1, 1)),
+ (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)),
+ (slice(i0, i1 + 1, 1), self.j1 + 1, slice(k0, k1, 1)), self.Hx],
+ ['hz_front_1', True,
+ (slice(i0, i1, 1), self.j0 - 2, slice(k0, k1 + 1, 1)),
+ (slice(i0, i1, 1), self.j0 - 1, slice(k0, k1 + 1, 1)),
+ (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)),
+ (slice(i0, i1, 1), self.j0 + 1, slice(k0, k1 + 1, 1)), self.Hz],
+ ['hz_back_1', True,
+ (slice(i0, i1, 1), self.j1 - 2, slice(k0, k1 + 1, 1)),
+ (slice(i0, i1, 1), self.j1 - 1, slice(k0, k1 + 1, 1)),
+ (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)),
+ (slice(i0, i1, 1), self.j1 + 1, slice(k0, k1 + 1, 1)), self.Hz],
+ ['hx_bottom_1', False,
+ # check these indexes
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 2),
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 1),
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0),
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 + 1), self.Hx],
+ ['hx_top_1', False,
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 2),
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 1),
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1),
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 + 1), self.Hx],
+ ['hy_bottom_1', True,
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 2),
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 1),
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0),
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 + 1), self.Hy],
+ ['hy_top_1', True,
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 2),
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 1),
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1),
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 + 1), self.Hy]
+ ]
+
+ for obj in slices:
+ sliced_field = obj[-1][obj[2]]
+ obj[1] = self.create_interpolated_coords(obj[1], sliced_field)
+
+ self.magnetic_slices = slices
+
+ def initialize_electric_slices_array(self):
+ # Extend the region sliced from the main grid by 1 cell.
+ # this allows more accurate interpolation for the outernodes
+ #i0 = self.i0 - 1
+ #j0 = self.j0 - 1
+ #k0 = self.k0 - 1
+ #i1 = self.i1 + 1
+ #j1 = self.j1 + 1
+ #k1 = self.k1 + 1
+
+ # not extended
+ i0 = self.i0
+ j0 = self.j0
+ k0 = self.k0
+ i1 = self.i1
+ j1 = self.j1
+ k1 = self.k1
+
+ # Spatially interpolate nodes
+ slices = [
+ ['ex_front_1', True,
+ (slice(i0, i1, 1), self.j0 - 1, slice(k0, k1 + 1, 1)),
+ (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)),
+ (slice(i0, i1, 1), self.j0 + 1, slice(k0, k1 + 1, 1)),
+ self.Ex],
+ ['ex_back_1', True,
+ (slice(i0, i1, 1), self.j1 - 1, slice(k0, k1 + 1, 1)),
+ (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)),
+ (slice(i0, i1, 1), self.j1 + 1, slice(k0, k1 + 1, 1)),
+ self.Ex],
+ ['ez_front_1', False,
+ (slice(i0, i1 + 1, 1), self.j0 - 1, slice(k0, k1, 1)),
+ (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)),
+ (slice(i0, i1 + 1, 1), self.j0 + 1, slice(k0, k1, 1)),
+ self.Ez],
+ ['ez_back_1', False,
+ (slice(i0, i1 + 1, 1), self.j1 - 1, slice(k0, k1, 1)),
+ (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)),
+ (slice(i0, i1 + 1, 1), self.j1 + 1, slice(k0, k1, 1)),
+ self.Ez],
+ ['ey_left_1', True,
+ (self.i0 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+ (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+ (self.i0 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+ self.Ey],
+ ['ey_right_1', True,
+ (self.i1 - 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+ (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+ (self.i1 + 1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
+ self.Ey],
+ ['ez_left_1', False,
+ (self.i0 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+ (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+ (self.i0 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+ self.Ez],
+ ['ez_right_1', False,
+ (self.i1 - 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+ (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+ (self.i1 + 1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
+ self.Ez],
+
+ ['ex_bottom_1', True,
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 - 1),
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0),
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0 + 1),
+ self.Ex],
+ ['ex_top_1', True,
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 - 1),
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1),
+ (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1 + 1),
+ self.Ex],
+ ['ey_bottom_1', False,
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 - 1),
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0),
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0 + 1),
+ self.Ey],
+ ['ey_top_1', False,
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 - 1),
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1),
+ (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1 + 1),
+ self.Ey]
+ ]
+
+ for obj in slices:
+ sliced_field = obj[-1][obj[2]]
+ obj[1] = self.create_interpolated_coords(obj[1], sliced_field)
+
+ self.electric_slices = slices
+
+ def create_interpolated_coords(self, mid, field):
+
+ n_x = field.shape[0]
+ n_y = field.shape[1]
+
+ if mid:
+ x = np.arange(0.5, n_x, 1.0)
+ z = np.arange(0, n_y, 1.0)
+
+ # Coordinates that require interpolated values
+ x_sg = np.linspace(self.d, n_x - self.d, n_x * self.ratio)
+ z_sg = np.linspace(0, n_y - 1, (n_y - 1) * self.ratio + 1)
+
+ else:
+ x = np.arange(0, n_x, 1.0)
+ z = np.arange(0.5, n_y, 1.0)
+
+ # Coordinates that require interpolated values
+ x_sg = np.linspace(0, n_x - 1, (n_x - 1) * self.ratio + 1)
+ z_sg = np.linspace(self.d, n_y - self.d, n_y * self.ratio)
+
+ return (x, z, x_sg, z_sg)
+
+ def update_magnetic(self):
+
+ # Copy previous time step magnetic field values to the previous
+ # time step variables
+ self.update_previous_timestep_fields(self.fn_m)
+
+ for obj in self.magnetic_slices:
+
+ # Grab the main grid fields used to interpolate across the IS
+ # f = self.Hi[slice]
+
+ f_u_1 = obj[-1][obj[2]]
+ f_u_2 = obj[-1][obj[3]]
+ f_u_3 = obj[-1][obj[4]]
+ f_u_4 = obj[-1][obj[5]]
+
+ f_1 = 0.25 * f_u_1 + 0.5 * f_u_2 + 0.25 * f_u_3
+ f_2 = 0.25 * f_u_2 + 0.5 * f_u_3 + 0.25 * f_u_4
+
+ #f_1 = obj[-1][obj[2]]
+ #f_2 = obj[-1][obj[3]]
+ if ('left' in obj[0] or
+ 'bottom' in obj[0] or
+ 'front' in obj[0]):
+ w = self.l_weight
+ else:
+ w = self.r_weight
+ c1, c2 = calculate_weighting_coefficients(w, self.ratio)
+ # transverse interpolated h field
+ f_t = c1 * f_1 + c2 * f_2
+
+ # interpolate over a fine grid
+ f_i = self.interpolate_to_sub_grid(f_t, obj[1])
+
+ if f_i == f_t:
+ raise ValueError
+
+ if obj[0] == 'hz_front_1' or obj[0] == 'hz_back_1':
+
+ pre_0 = int(f_i.shape[0] / 2 - 1 - self.ratio // 2)
+ pre_1 = int(f_i.shape[0] / 2 + self.ratio // 2)
+
+ f_i[pre_0 + 1:pre_0 + 1 + self.ratio // 2, :] = np.stack((f_i[pre_0, :], f_i[pre_0, :]))
+ f_i[pre_1 - self.ratio // 2:pre_1, :] = np.stack((f_i[pre_1, :], f_i[pre_1, :]))
+
+
+ # discard the outer nodes only required for interpolation
+ #f = f_i[self.ratio:-self.ratio, self.ratio:-self.ratio]
+ f = f_i
+ setattr(self, obj[0], f)
+
+ def update_electric(self):
+
+ self.update_previous_timestep_fields(self.fn_e)
+
+ for obj in self.electric_slices:
+
+ f_u_1 = obj[-1][obj[2]]
+ f_u_2 = obj[-1][obj[3]]
+ f_u_3 = obj[-1][obj[4]]
+
+ f_m = 0.25 * f_u_1 + 0.5 * f_u_2 + 0.25 * f_u_3
+
+ f_i = self.interpolate_to_sub_grid(f_m, obj[1])
+
+ # correction interpolated points at the pec traversal point
+ if obj[0] == 'ex_front_1' or obj[0] == 'ex_back_1':
+
+ pre_0 = int(f_i.shape[0] / 2 - 1 - self.ratio // 2)
+ pre_1 = int(f_i.shape[0] / 2 + self.ratio // 2)
+
+ f_i[pre_0 + 1:pre_0 + 1 + self.ratio // 2, :] = np.stack((f_i[pre_0, :], f_i[pre_0, :]))
+ f_i[pre_1 - self.ratio // 2:pre_1, :] = np.stack((f_i[pre_1, :], f_i[pre_1, :]))
+
+ f = f_i
+ #f = f_i[self.ratio:-self.ratio, self.ratio:-self.ratio]
+ setattr(self, obj[0], f)
+
+ def interpolate_to_sub_grid(self, field, coords):
+ x, z, x_sg, z_sg = coords
+ ex_t = np.transpose(field)
+ f = interpolate.interp2d(x, z, ex_t, kind=self.interpolation)
+ #f = interpolate.RectBivariateSpline(x, z, field)
+ ex_inter_t = f(x_sg, z_sg)
+ ex_inter = np.transpose(ex_inter_t)
+ #ex_inter = ex_inter_t
+
+ return ex_inter
+
+
+class PlaneError(Exception):
+ pass
diff --git a/gprMax/subgrids/solver.py b/gprMax/subgrids/solver.py
new file mode 100644
index 00000000..4c164d25
--- /dev/null
+++ b/gprMax/subgrids/solver.py
@@ -0,0 +1,355 @@
+from ..cython.fields_updates_normal import update_electric
+from ..cython.fields_updates_normal import update_magnetic
+from ..fields_outputs import store_outputs
+from ..utilities import get_terminal_width
+from ..exceptions import GeneralError
+
+from .subgrid_hsg import SubGridHSG
+from .precursor_nodes import PrecursorNodes as PrecursorNodesHSG
+from .precursor_nodes_filtered import PrecursorNodes as PrecursorNodesFilteredHSG
+
+from tqdm import tqdm
+from time import perf_counter
+
+import os
+import sys
+
+
+def create_solver(G):
+ """Return the solver for the given subgrids."""
+ updaters = []
+
+ for sg in G.subgrids:
+ print(sg)
+ sg_type = type(sg)
+ if sg_type == SubGridHSG and sg.filter:
+ precursors = PrecursorNodesFilteredHSG(G, sg)
+ elif sg_type == SubGridHSG and not sg.filter:
+ precursors = PrecursorNodesHSG(G, sg)
+ else:
+ raise GeneralError(str(sg) + ' is not a subgrid type')
+
+ sgu = SubgridUpdater(sg, precursors, G)
+ updaters.append(sgu)
+
+ solver = SubGridSolver(G, updaters)
+ return solver
+
+
+class SubGridSolver:
+ """Solver for subgridding simulations."""
+
+ """Class to call the various update methods required for an HSG-Subgrid simulation.
+ Multiple subgrids can be updated by adding more subgrid_updater objects to the subgrid_updater
+ array.
+ """
+
+ def __init__(self, G, subgrid_updaters, hsg=True):
+ """
+ Args:
+ G (G): Grid class instance - holds essential parameters
+ describing the model.
+ subgrid_updaters: (list): list of subgrid_updaters used for updating
+ the subgrids
+ iterations (int): number of iterations for the simulation.
+ hsg (bool): HSG methods for subgrids will not be called if False.
+
+ """
+
+ self.G = G
+ self.grids = [self.G] + G.subgrids
+ self.subgrid_updaters = subgrid_updaters
+ self.materials = G.materials
+ self.abs_time = 0
+ self.hsg = hsg
+
+ def store_outputs(self):
+ """Method to store field outputs for all grid for each main grid iteration."""
+ for grid in self.grids:
+ store_outputs(self.G.iteration, grid.Ex, grid.Ey, grid.Ez,
+ grid.Hx, grid.Hy, grid.Hz, grid)
+
+ def store_snapshots(self):
+ """Store any snapshots."""
+ for snap in self.G.snapshots:
+ if snap.time == self.G.iteration + 1:
+ snap.store(self.G)
+
+ 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 solve(self, iterations):
+ """Run timestepping."""
+ tsolvestart = perf_counter()
+ self.iterations = iterations
+ # for time step in range(self.G.iterations):
+
+ # The main grid FDTD loop
+ for iteration in self.iterations:
+
+ # Keep track of the index. Required for saving output correctly
+ self.G.iteration = iteration
+
+ # Write any snapshots of the E, H, I fields/currents
+ self.write_snapshots(iteration)
+
+ self.store_outputs()
+ # Update main grid electric field components including sources
+ self.update_magnetic()
+ # Update the fields in the subgrids / main grid
+ if self.hsg:
+ self.hsg_2()
+
+ # Update main grid electric field components including sources
+ self.update_electric()
+ # Update the fields in the subgrids / main grid
+ if self.hsg:
+ self.hsg_1()
+
+ # Return the elapsed time
+ tsolve = perf_counter() - tsolvestart
+
+ return tsolve
+
+ def update_electric(self):
+ """Method to update E fields, PML and electric sources."""
+ # All materials are non-dispersive so do standard update
+ G = self.G
+
+ update_electric(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsE, G.ID,
+ G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
+
+ # Update electric field components with the PML correction
+ for pml in G.pmls:
+ pml.update_electric(G)
+ # Update electric field components from sources (update any Hertzian dipole sources last)
+ for source in G.voltagesources + G.transmissionlines + G.hertziandipoles:
+ source.update_electric(G.iteration, G.updatecoeffsE, G.ID, G.Ex,
+ G.Ey, G.Ez, G)
+
+ def update_magnetic(self):
+ """Method to update H fields, PML and magnetic sources."""
+ # Update magnetic field components
+ G = self.G
+
+ update_magnetic(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsH, G.ID,
+ G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
+
+ # Update magnetic field components with the PML correction
+ for pml in G.pmls:
+ pml.update_magnetic(G)
+
+ # Update magnetic field components from sources
+ for source in G.transmissionlines + G.magneticdipoles:
+ source.update_magnetic(G.iteration, G.updatecoeffsH, G.ID,
+ G.Hx, G.Hy, G.Hz, G)
+
+ def write_snapshots(self, iteration):
+ # Write any snapshots to file
+ for i, snap in enumerate(self.G.snapshots):
+ if snap.time == iteration + 1:
+ snapiters = 36 * (((snap.xf - snap.xs) / snap.dx) * ((snap.yf - snap.ys) / snap.dy) * ((snap.zf - snap.zs) / snap.dz))
+ pbar = tqdm(total=snapiters, leave=False, unit='byte', unit_scale=True, desc=' Writing snapshot file {} of {}, {}'.format(i + 1, len(self.G.snapshots), os.path.split(snap.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=self.G.tqdmdisable)
+
+ # Use this call to print out main grid and subgrids
+ snap.write_vtk_imagedata(self.G.Ex, self.G.Ey, self.G.Ez, self.G.Hx, self.G.Hy, self.G.Hz, self.G, pbar, sub_grids=self.G.subgrids)
+
+ # Use this call to print out the standard grid without subgrid
+ # snap.write_vtk_imagedata(self.G.Ex, self.G.Ey, self.G.Ez, self.G.Hx, self.G.Hy, self.G.Hz, self.G, pbar)
+
+ # Use this call to print out only the subgrid - use in combination with commented code in .multi_cmds/snapshots.py
+ # snap.write_vtk_imagedata_fast(self.grid)
+ pbar.close()
+
+
+class SubgridUpdater:
+ """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.
+ """
+ self.subgrid = subgrid
+ self.precursors = precursors
+ self.G = G
+ self.source_iteration = 0
+
+ def hsg_1(self):
+ """This is the first half of the subgrid update. Takes the time step
+ up to the main grid magnetic update"""
+ G = self.G
+ sub_grid = self.subgrid
+ precursors = self.precursors
+
+ precursors.update_electric()
+
+ 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
+ update_electric(sub_grid.nx,
+ sub_grid.ny,
+ sub_grid.nz,
+ G.nthreads,
+ sub_grid.updatecoeffsE,
+ sub_grid.ID,
+ sub_grid.Ex,
+ sub_grid.Ey,
+ sub_grid.Ez,
+ sub_grid.Hx,
+ sub_grid.Hy,
+ sub_grid.Hz)
+ for pml in sub_grid.pmls:
+ pml.update_electric(sub_grid)
+ precursors.interpolate_magnetic_in_time(int(m + sub_grid.ratio / 2 - 0.5))
+ sub_grid.update_electric_is(precursors)
+
+ self.update_sub_grid_electric_sources()
+
+ # STD update, interpolate inc. field in time, apply correction
+ update_magnetic(sub_grid.nx,
+ sub_grid.ny,
+ sub_grid.nz,
+ G.nthreads,
+ sub_grid.updatecoeffsH,
+ sub_grid.ID,
+ sub_grid.Ex,
+ sub_grid.Ey,
+ sub_grid.Ez,
+ sub_grid.Hx,
+ sub_grid.Hy,
+ sub_grid.Hz)
+ 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)
+ update_electric(sub_grid.nx,
+ sub_grid.ny,
+ sub_grid.nz,
+ G.nthreads,
+ sub_grid.updatecoeffsE,
+ sub_grid.ID,
+ sub_grid.Ex,
+ sub_grid.Ey,
+ sub_grid.Ez,
+ sub_grid.Hx,
+ sub_grid.Hy,
+ sub_grid.Hz)
+ for pml in sub_grid.pmls:
+ pml.update_electric(sub_grid)
+ precursors.calc_exact_magnetic_in_time()
+ sub_grid.update_electric_is(precursors)
+ self.update_sub_grid_electric_sources()
+
+ sub_grid.update_electric_os(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"""
+ G = self.G
+ sub_grid = self.subgrid
+ precursors = self.precursors
+
+ precursors.update_magnetic()
+
+ upper_m = int(sub_grid.ratio / 2 - 0.5)
+
+ for m in range(1, upper_m + 1):
+
+ update_magnetic(sub_grid.nx,
+ sub_grid.ny,
+ sub_grid.nz,
+ G.nthreads,
+ sub_grid.updatecoeffsH,
+ sub_grid.ID,
+ sub_grid.Ex,
+ sub_grid.Ey,
+ sub_grid.Ez,
+ sub_grid.Hx,
+ sub_grid.Hy,
+ sub_grid.Hz)
+
+ 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)
+ update_electric(sub_grid.nx,
+ sub_grid.ny,
+ sub_grid.nz,
+ G.nthreads,
+ sub_grid.updatecoeffsE,
+ sub_grid.ID,
+ sub_grid.Ex,
+ sub_grid.Ey,
+ sub_grid.Ez,
+ sub_grid.Hx,
+ sub_grid.Hy,
+ sub_grid.Hz)
+
+ for pml in sub_grid.pmls:
+ pml.update_electric(sub_grid)
+
+ precursors.interpolate_magnetic_in_time(m)
+ sub_grid.update_electric_is(precursors)
+ self.update_sub_grid_electric_sources()
+
+ update_magnetic(sub_grid.nx,
+ sub_grid.ny,
+ sub_grid.nz,
+ G.nthreads,
+ sub_grid.updatecoeffsH,
+ sub_grid.ID,
+ sub_grid.Ex,
+ sub_grid.Ey,
+ sub_grid.Ez,
+ sub_grid.Hx,
+ sub_grid.Hy,
+ sub_grid.Hz)
+ 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(G)
+
+ def update_sub_grid_electric_sources(self):
+ """Update any electric sources in the subgrid"""
+ sg = self.subgrid
+ 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.subgrid
+ 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/hsg.py b/gprMax/subgrids/subgrid_hsg.py
similarity index 58%
rename from gprMax/subgrids/hsg.py
rename to gprMax/subgrids/subgrid_hsg.py
index 3327b40a..8aa1d1a1 100644
--- a/gprMax/subgrids/hsg.py
+++ b/gprMax/subgrids/subgrid_hsg.py
@@ -1,25 +1,7 @@
-# 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 .base import SubGridBase
from ..cython.fields_updates_hsg import cython_update_is
-from ..cython.fields_updates_hsg import cython_update_electric_os
from ..cython.fields_updates_hsg import cython_update_magnetic_os
+from ..cython.fields_updates_hsg import cython_update_electric_os
from ..utilities import human_size
from colorama import init, Fore, Style
@@ -34,89 +16,98 @@ class SubGridHSG(SubGridBase):
super().__init__(**kwargs)
self.gridtype = SubGridHSG.gridtype
- def memory_usage(self):
- pass
-
def update_magnetic_is(self, precursors):
- """Apply the incident field corrections to the subgrid fields at the IS.
-
- Args:
- precursors (PrecursorNodes): PrecursorNodes are the incident
- fields calculated at the correct time and position.
+ """Update the subgrid nodes at the IS with the currents derived
+ from the main grid.
+ Args: nwl, nwm, nwn, face, field, inc_field, lookup_id, sign, mod, co
"""
- # bottom and top faces
+ # Hz = c0Hz - c1Ey + c2Ex
+ # Hy = c0Hy - c3Ex + c1Ez
+ # Hx = c0Hx - c2Ez + c3Ey
+ # bottom and top
cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx, self.nwy + 1, self.nwz, 1, self.Hy, precursors.ex_bottom, precursors.ex_top, self.IDlookup['Hy'], 1, -1, 3, self.nthreads)
cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwy, self.nwz, 1, self.Hx, precursors.ey_bottom, precursors.ey_top, self.IDlookup['Hx'], -1, 1, 3, self.nthreads)
- # left and right faces
+ # left and right
cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwy, self.nwz + 1, self.nwx, 2, self.Hz, precursors.ey_left, precursors.ey_right, self.IDlookup['Hz'], 1, -1, 1, self.nthreads)
cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwy + 1, self.nwz, self.nwx, 2, self.Hy, precursors.ez_left, precursors.ez_right, self.IDlookup['Hy'], -1, 1, 1, self.nthreads)
- # front and back faces
+ # front and back
cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx, self.nwz + 1, self.nwy, 3, self.Hz, precursors.ex_front, precursors.ex_back, self.IDlookup['Hz'], -1, 1, 2, self.nthreads)
cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwz, self.nwy, 3, self.Hx, precursors.ez_front, precursors.ez_back, self.IDlookup['Hx'], 1, -1, 2, self.nthreads)
def update_electric_is(self, precursors):
- """Apply the incident field corrections to the subgrid fields at the IS.
+ # Args: nwl, nwm, nwn, face, field, inc_field, lookup_id, sign, mod, co
- Args:
- precursors (PrecursorNodes): PrecursorNodes are the incident
- fields calculated at the correct time and position.
- """
- # bottom and top faces
+ # Ex = c0(Ex) + c2(dHz) - c3(dHy)
+ # Ey = c0(Ey) + c3(dHx) - c1(dHz)
+ # Ez = c0(Ez) + c1(dHy) - c2(dHx)
+
+ # bottom and top
cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx, self.nwy + 1, self.nwz, 1, self.Ex, precursors.hy_bottom, precursors.hy_top, self.IDlookup['Ex'], 1, -1, 3, self.nthreads)
cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx + 1, self.nwy, self.nwz, 1, self.Ey, precursors.hx_bottom, precursors.hx_top, self.IDlookup['Ey'], -1, 1, 3, self.nthreads)
- # left and right faces
+ # left and right
cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwy, self.nwz + 1, self.nwx, 2, self.Ey, precursors.hz_left, precursors.hz_right, self.IDlookup['Ey'], 1, -1, 1, self.nthreads)
cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwy + 1, self.nwz, self.nwx, 2, self.Ez, precursors.hy_left, precursors.hy_right, self.IDlookup['Ez'], -1, 1, 1, self.nthreads)
- # front and back faces
+ # front and back
cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx, self.nwz + 1, self.nwy, 3, self.Ex, precursors.hz_front, precursors.hz_back, self.IDlookup['Ex'], -1, 1, 2, self.nthreads)
cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx + 1, self.nwz, self.nwy, 3, self.Ez, precursors.hx_front, precursors.hx_back, self.IDlookup['Ez'], 1, -1, 2, self.nthreads)
def update_electric_os(self, main_grid):
- """Apply the incident field corrections to the main grid fields at the OS.
+ i_l = self.i0 - self.is_os_sep
+ i_u = self.i1 + self.is_os_sep
+ j_l = self.j0 - self.is_os_sep
+ j_u = self.j1 + self.is_os_sep
+ k_l = self.k0 - self.is_os_sep
+ k_u = self.k1 + self.is_os_sep
- Main grid fields are collocated with subgrid fields therefore
- precursor fields are not required.
+ # Args: sub_grid, normal, l_l, l_u, m_l, m_u, n_l, n_u, nwn, lookup_id, field, inc_field, co, sign_n, sign_f
- Args:
- main_grid (FDTDGrid): main fdtd solver grid
- """
- # front and back
- cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, self.i_l, self.i_u, self.k_l, self.k_u + 1, self.j_l, self.j_u, self.nwy, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hz, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
- cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, self.i_l, self.i_u + 1, self.k_l, self.k_u, self.j_l, self.j_u, self.nwy, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hx, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+ # Form of FDTD update equations for E
+ # Ex = c0(Ex) + c2(dHz) - c3(dHy)
+ # Ey = c0(Ey) + c3(dHx) - c1(dHz)
+ # Ez = c0(Ez) + c1(dHy) - c2(dHx)
+
+ cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, i_l, i_u, k_l, k_u + 1, j_l, j_u, self.nwy, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hz, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+ cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, i_l, i_u + 1, k_l, k_u, j_l, j_u, self.nwy, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hx, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
# Left and Right
- cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, self.j_l, self.j_u, self.k_l, self.k_u + 1, self.i_l, self.i_u, self.nwx, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hz, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
- cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, self.j_l, self.j_u + 1, self.k_l, self.k_u, self.i_l, self.i_u, self.nwx, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hy, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+ cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, j_l, j_u, k_l, k_u + 1, i_l, i_u, self.nwx, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hz, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+ cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, j_l, j_u + 1, k_l, k_u, i_l, i_u, self.nwx, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hy, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
# Bottom and Top
- cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, self.i_l, self.i_u, self.j_l, self.j_u + 1, self.k_l, self.k_u, self.nwz, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hy, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
- cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, self.i_l, self.i_u + 1, self.j_l, self.j_u, self.k_l, self.k_u, self.nwz, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hx, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+ cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, i_l, i_u, j_l, j_u + 1, k_l, k_u, self.nwz, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hy, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+ cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l, k_u, self.nwz, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hx, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
def update_magnetic_os(self, main_grid):
- """Apply the incident field corrections to the main grid fields at the OS.
- Main grid fields are collocated with subgrid fields therefore
- precursor fields are not required.
+ i_l = self.i0 - self.is_os_sep
+ i_u = self.i1 + self.is_os_sep
+ j_l = self.j0 - self.is_os_sep
+ j_u = self.j1 + self.is_os_sep
+ k_l = self.k0 - self.is_os_sep
+ k_u = self.k1 + self.is_os_sep
- Args:
- main_grid (FDTDGrid): main fdtd solver grid
- """
+ # Form of FDTD update equations for H
+ # Hz = c0Hz - c1Ey + c2Ex
+ # Hy = c0Hy - c3Ex + c1Ez
+ # Hx = c0Hx - c2Ez + c3Ey
+
+ # Args: sub_grid, normal, l_l, l_u, m_l, m_u, n_l, n_u, nwn, lookup_id, field, inc_field, co, sign_n, sign_f):
# Front and back
- cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, self.i_l, self.i_u, self.k_l, self.k_u + 1, self.j_l - 1, self.j_u, self.nwy, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ex, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
- cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, self.i_l, self.i_u + 1, self.k_l, self.k_u, self.j_l - 1, self.j_u, self.nwy, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ez, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+ cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, i_l, i_u, k_l, k_u + 1, j_l - 1, j_u, self.nwy, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ex, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+ cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, i_l, i_u + 1, k_l, k_u, j_l - 1, j_u, self.nwy, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ez, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
# Left and Right
- cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, self.j_l, self.j_u, self.k_l, self.k_u + 1, self.i_l - 1, self.i_u, self.nwx, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ey, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
- cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, self.j_l, self.j_u + 1, self.k_l, self.k_u, self.i_l - 1, self.i_u, self.nwx, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ez, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+ cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, j_l, j_u, k_l, k_u + 1, i_l - 1, i_u, self.nwx, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ey, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+ cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, j_l, j_u + 1, k_l, k_u, i_l - 1, i_u, self.nwx, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ez, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
# bottom and top
- cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, self.i_l, self.i_u, self.j_l, self.j_u + 1, self.k_l - 1, self.k_u, self.nwz, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ex, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
- 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)
+ cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u, j_l, j_u + 1, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ex, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, main_grid.nthreads)
+ cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l - 1, 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):
diff --git a/gprMax/subgrids/updates.py b/gprMax/subgrids/updates.py
index 2fb957f4..be3759b7 100644
--- a/gprMax/subgrids/updates.py
+++ b/gprMax/subgrids/updates.py
@@ -16,6 +16,9 @@
# You should have received a copy of the GNU General Public License
# along with gprMax. If not, see .
from ..updates import CPUUpdates
+from gprMax.utilities import timer
+from ..cython.fields_updates_normal import update_electric
+from ..cython.fields_updates_normal import update_magnetic
def create_subgrid_updates(G):
@@ -41,6 +44,78 @@ def create_subgrid_updates(G):
return updates
+class SubgridSolver:
+
+ """Generic solver for Update objects"""
+
+ def __init__(self, updates):
+ """Context for the model to run in. Sub-class this with contexts
+ i.e. an MPI context.
+
+ Args:
+ updates (Updates): updates contains methods to run FDTD algorithm
+ iterator (iterator): can be range() or tqdm()
+ """
+ self.updates = updates
+
+ def get_G(self):
+ return self.updates.grid
+
+ def update_electric_1(self):
+ """Method to update E fields, PML and electric sources"""
+ G = self.get_G()
+ update_electric(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
+
+ # Update electric field components with the PML correction
+ for pml in G.pmls:
+ pml.update_electric(G)
+ # Update electric field components from sources (update any Hertzian dipole sources last)
+ for source in G.voltagesources + G.transmissionlines + G.hertziandipoles:
+ source.update_electric(G.iteration, G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G)
+
+ def update_electric_2(self):
+ return
+
+ def update_magnetic(self):
+ """Method to update H fields, PML and magnetic sources"""
+ # Update magnetic field components
+ G = self.get_G()
+
+ update_magnetic(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsH, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
+
+ # Update magnetic field components with the PML correction
+ for pml in G.pmls:
+ pml.update_magnetic(G)
+
+ # Update magnetic field components from sources
+ for source in G.transmissionlines + G.magneticdipoles:
+ source.update_magnetic(G.iteration, G.updatecoeffsH, G.ID, G.Hx, G.Hy, G.Hz, G)
+
+ def solve(self, iterator):
+ """Time step the FDTD model."""
+ tsolvestart = timer()
+
+ for iteration in iterator:
+ # Update main grid electric field components including sources
+ self.updates.store_outputs(iteration)
+ self.update_electric_1()
+ #self.updates.update_electric_pml()
+ #self.updates.update_electric_sources(iteration)
+ # Update the fields in the subgrids / main grid
+ self.updates.hsg_1()
+ # apply 2nd dispersive update after OS updates
+ #self.updates.update_electric_b()
+ # Update main grid electric field components including sources
+ self.update_magnetic()
+ #self.updates.update_magnetic_pml()
+ #self.updates.update_magnetic_sources(iteration)
+ # Update the fields in the subgrids / main grid
+ self.updates.hsg_2()
+
+ tsolve = timer() - tsolvestart
+ return tsolve
+
+
class SubgridUpdates(CPUUpdates):
"""Top level subgrid updater. Manages the collection of subgrids."""
@@ -48,6 +123,11 @@ class SubgridUpdates(CPUUpdates):
super().__init__(G)
self.subgrid_updaters = subgrid_updaters
+ def store_outputs(self, iteration):
+ super().store_outputs(iteration)
+ for updater in self.subgrid_updaters:
+ updater.store_outputs(iteration)
+
def hsg_1(self):
"""Method to update the subgrids over the first phase"""
for sg_updater in self.subgrid_updaters:
@@ -58,14 +138,6 @@ class SubgridUpdates(CPUUpdates):
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)
@@ -94,33 +166,51 @@ class SubgridUpdater(CPUUpdates):
self.source_iteration = 0
self.G = G
+ def update_electric_2(self):
+ return
+
+ def update_electric_1(self):
+ subgrid = self.grid
+ # Update electric field components
+ # All materials are non-dispersive so do standard update
+ update_electric(subgrid.nx, subgrid.ny, subgrid.nz, subgrid.nthreads, subgrid.updatecoeffsE, subgrid.ID, subgrid.Ex, subgrid.Ey, subgrid.Ez, subgrid.Hx, subgrid.Hy, subgrid.Hz)
+
def hsg_1(self):
- """This is the first half of the subgrid update. Takes the time step
- up to the main grid magnetic update"""
+ G = self.G
+ sub_grid = self.grid
+
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()
+ self.update_electric_1()
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()
+ self.update_electric_2()
# STD update, interpolate inc. field in time, apply correction
- self.update_magnetic()
+ update_magnetic(sub_grid.nx,
+ sub_grid.ny,
+ sub_grid.nz,
+ G.nthreads,
+ sub_grid.updatecoeffsH,
+ sub_grid.ID,
+ sub_grid.Ex,
+ sub_grid.Ey,
+ sub_grid.Ez,
+ sub_grid.Hx,
+ sub_grid.Hy,
+ sub_grid.Hz)
for pml in sub_grid.pmls:
pml.update_magnetic(sub_grid)
@@ -129,7 +219,7 @@ class SubgridUpdater(CPUUpdates):
self.update_sub_grid_magnetic_sources()
# store_outputs(self.grid)
- self.update_electric_a()
+ self.update_electric_1()
for pml in sub_grid.pmls:
pml.update_electric(sub_grid)
@@ -137,13 +227,12 @@ class SubgridUpdater(CPUUpdates):
precursors.calc_exact_magnetic_in_time()
sub_grid.update_electric_is(precursors)
- self.update_electric_b()
+ self.update_electric_2()
- sub_grid.update_electric_os(self.G)
+ sub_grid.update_electric_os(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"""
+ G = self.G
sub_grid = self.grid
precursors = self.precursors
@@ -153,7 +242,18 @@ class SubgridUpdater(CPUUpdates):
for m in range(1, upper_m + 1):
- self.update_magnetic()
+ update_magnetic(sub_grid.nx,
+ sub_grid.ny,
+ sub_grid.nz,
+ G.nthreads,
+ sub_grid.updatecoeffsH,
+ sub_grid.ID,
+ sub_grid.Ex,
+ sub_grid.Ey,
+ sub_grid.Ez,
+ sub_grid.Hx,
+ sub_grid.Hy,
+ sub_grid.Hz)
for pml in sub_grid.pmls:
pml.update_magnetic(sub_grid)
@@ -163,7 +263,7 @@ class SubgridUpdater(CPUUpdates):
self.update_sub_grid_magnetic_sources()
# store_outputs(self.grid)
- self.update_electric_a()
+ self.update_electric_1()
for pml in sub_grid.pmls:
pml.update_electric(sub_grid)
@@ -171,9 +271,22 @@ class SubgridUpdater(CPUUpdates):
precursors.interpolate_magnetic_in_time(m)
sub_grid.update_electric_is(precursors)
- self.update_electric_b()
+ self.update_electric_2()
- self.update_magnetic()
+ #sub_grid.update_magnetic_os(G)
+
+ update_magnetic(sub_grid.nx,
+ sub_grid.ny,
+ sub_grid.nz,
+ G.nthreads,
+ sub_grid.updatecoeffsH,
+ sub_grid.ID,
+ sub_grid.Ex,
+ sub_grid.Ey,
+ sub_grid.Ez,
+ sub_grid.Hx,
+ sub_grid.Hy,
+ sub_grid.Hz)
for pml in sub_grid.pmls:
pml.update_magnetic(sub_grid)
@@ -182,7 +295,7 @@ class SubgridUpdater(CPUUpdates):
sub_grid.update_magnetic_is(precursors)
self.update_sub_grid_magnetic_sources()
- sub_grid.update_magnetic_os(self.G)
+ sub_grid.update_magnetic_os(G)
def update_sub_grid_electric_sources(self):
"""Update any electric sources in the subgrid"""
diff --git a/gprMax/subgrids/user_objects.py b/gprMax/subgrids/user_objects.py
index 624115c2..34104592 100644
--- a/gprMax/subgrids/user_objects.py
+++ b/gprMax/subgrids/user_objects.py
@@ -1,11 +1,10 @@
-from .hsg import SubGridHSG as SubGridHSGUser
+from .subgrid_hsg import SubGridHSG as SubGridHSGUser
from .multi import ReferenceRx as ReferenceRxUser
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 gprMax import config
from copy import copy
@@ -39,23 +38,33 @@ class SubGridBase(UserObjectMulti):
def set_discretisation(self, sg, grid):
"""Set the spatial discretisation."""
- sg.dl = grid.dl / sg.ratio
- sg.dx, sg.dy, sg.dz = grid.dl
+ sg.dx = grid.dx / sg.ratio
+ sg.dy = grid.dy / sg.ratio
+ sg.dz = grid.dz / sg.ratio
+ sg.dl = np.array([sg.dx, sg.dy, sg.dz])
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
- # IS indices
- sg.i0, sg.j0, sg.k0 = p1
- sg.i1, sg.j1, sg.k1 = 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)
- # 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
+ # 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)
+
+ 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)
+
+ def set_name(self):
+ self.name = self.kwargs['id']
- # 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."""
@@ -73,9 +82,6 @@ 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']
@@ -83,8 +89,6 @@ 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)
@@ -98,9 +102,20 @@ 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)
+ self.set_name()
# Copy a reference for the main grid to the sub grid
sg.parent_grid = grid
@@ -133,9 +148,11 @@ class SubGridHSG(SubGridBase):
ratio=3,
ID='',
is_os_sep=3,
+ pml_separation=4,
subgrid_pml_thickness=6,
- interpolation=3,
+ interpolation='linear',
loss_mechanism=False,
+ loss_factor=False,
filter=True,
**kwargs):
"""Constructor."""
@@ -146,10 +163,12 @@ class SubGridHSG(SubGridBase):
kwargs['ratio'] = ratio
kwargs['ID'] = ID
kwargs['is_os_sep'] = is_os_sep
- kwargs['pml_separation'] = ratio // 2 + 2
+ kwargs['pml_separation'] = pml_separation
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
@@ -158,10 +177,12 @@ class SubGridHSG(SubGridBase):
def create(self, grid, uip):
sg = SubGridHSGUser(**self.kwargs)
self.setup(sg, grid, uip)
- if config.general['messages']:
+ if config.is_messages():
print(sg)
return sg
+
+
class ReferenceRx(Rx):
"""ReferenceRx User Object."""
diff --git a/gprMax/utilities.py b/gprMax/utilities.py
index ef92238e..945a67de 100644
--- a/gprMax/utilities.py
+++ b/gprMax/utilities.py
@@ -38,6 +38,7 @@ import numpy as np
from .exceptions import GeneralError
+
def get_terminal_width():
"""Get/set width of terminal being used.
@@ -421,8 +422,8 @@ def timer():
class Printer():
- def __init__(self, sim_config):
- self.printing = sim_config.general['messages']
+ def __init__(self, config):
+ self.printing = config.is_messages()
def print(self, str):
if self.printing: