你已经派生过 gprMax
镜像自地址
https://gitee.com/sunhf/gprMax.git
已同步 2025-08-07 23:14:03 +08:00
Whitespace cleanups.
这个提交包含在:
@@ -16,10 +16,12 @@
|
||||
# You should have received a copy of the GNU General Public License
|
||||
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
|
||||
class GeneralError(ValueError):
|
||||
"""Handles general errors. Subclasses the ValueError class."""
|
||||
pass
|
||||
|
||||
|
||||
class CmdInputError(ValueError):
|
||||
"""Handles errors in user specified commands. Subclasses the ValueError class."""
|
||||
pass
|
||||
|
@@ -76,7 +76,7 @@ class FractalSurface(object):
|
||||
self.fractalsurface = np.zeros(surfacedims, dtype=complextype)
|
||||
|
||||
# Positional vector at centre of array, scaled by weighting
|
||||
v1 = np.array([self.weighting[0]*(surfacedims[0])/2, self.weighting[1]*(surfacedims[1])/2])
|
||||
v1 = np.array([self.weighting[0] * (surfacedims[0]) / 2, self.weighting[1] * (surfacedims[1]) / 2])
|
||||
|
||||
# 2D array of random numbers to be convolved with the fractal function
|
||||
R = np.random.RandomState(self.seed)
|
||||
@@ -88,13 +88,13 @@ class FractalSurface(object):
|
||||
for i in range(surfacedims[0]):
|
||||
for j in range(surfacedims[1]):
|
||||
# Positional vector for current position
|
||||
v2 = np.array([self.weighting[0]*i, self.weighting[1]*j])
|
||||
v2 = np.array([self.weighting[0] * i, self.weighting[1] * j])
|
||||
rr = np.linalg.norm(v2 - v1)
|
||||
try:
|
||||
self.fractalsurface[i, j] = A[i, j] * 1/(rr**self.b)
|
||||
self.fractalsurface[i, j] = A[i, j] * 1 / (rr**self.b)
|
||||
except FloatingPointError:
|
||||
rr = 0.9
|
||||
self.fractalsurface[i, j] = A[i, j] * 1/(rr**self.b)
|
||||
self.fractalsurface[i, j] = A[i, j] * 1 / (rr**self.b)
|
||||
|
||||
# Shift the zero frequency component to the centre of the spectrum
|
||||
self.fractalsurface = np.fft.ifftshift(self.fractalsurface)
|
||||
@@ -104,7 +104,7 @@ class FractalSurface(object):
|
||||
fractalmin = np.amin(self.fractalsurface)
|
||||
fractalmax = np.amax(self.fractalsurface)
|
||||
fractalrange = fractalmax - fractalmin
|
||||
self.fractalsurface = self.fractalsurface * ((self.fractalrange[1] - self.fractalrange[0])/fractalrange) + self.fractalrange[0] - ((self.fractalrange[1] - self.fractalrange[0])/fractalrange) * fractalmin
|
||||
self.fractalsurface = self.fractalsurface * ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) + self.fractalrange[0] - ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * fractalmin
|
||||
|
||||
|
||||
class FractalVolume(object):
|
||||
@@ -118,7 +118,7 @@ class FractalVolume(object):
|
||||
"""
|
||||
|
||||
self.ID = None
|
||||
self.operatingonID= None
|
||||
self.operatingonID = None
|
||||
self.xs = xs
|
||||
self.xf = xf
|
||||
self.ys = ys
|
||||
@@ -146,7 +146,7 @@ class FractalVolume(object):
|
||||
self.fractalvolume = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=complextype)
|
||||
|
||||
# Positional vector at centre of array, scaled by weighting
|
||||
v1 = np.array([self.weighting[0]*(self.nx + 1)/2, self.weighting[1]*(self.ny + 1)/2, self.weighting[2]*(self.nz + 1)/2])
|
||||
v1 = np.array([self.weighting[0] * (self.nx + 1) / 2, self.weighting[1] * (self.ny + 1) / 2, self.weighting[2] * (self.nz + 1) / 2])
|
||||
|
||||
# 3D array of random numbers to be convolved with the fractal function
|
||||
R = np.random.RandomState(self.seed)
|
||||
@@ -159,13 +159,13 @@ class FractalVolume(object):
|
||||
for j in range(self.ny + 1):
|
||||
for k in range(self.nz + 1):
|
||||
# Positional vector for current position
|
||||
v2 = np.array([self.weighting[0]*i, self.weighting[1]*j, self.weighting[2]*k])
|
||||
v2 = np.array([self.weighting[0] * i, self.weighting[1] * j, self.weighting[2] * k])
|
||||
rr = np.linalg.norm(v2 - v1)
|
||||
try:
|
||||
self.fractalvolume[i, j, k] = A[i, j, k] * 1/(rr**self.b)
|
||||
self.fractalvolume[i, j, k] = A[i, j, k] * 1 / (rr**self.b)
|
||||
except FloatingPointError:
|
||||
rr = 0.9
|
||||
self.fractalvolume[i, j, k] = A[i, j, k] * 1/(rr**self.b)
|
||||
self.fractalvolume[i, j, k] = A[i, j, k] * 1 / (rr**self.b)
|
||||
|
||||
# Shift the zero frequency component to the centre of the spectrum
|
||||
self.fractalvolume = np.fft.ifftshift(self.fractalvolume)
|
||||
|
@@ -153,7 +153,7 @@ class GeometryView(object):
|
||||
vtk_numpoints = (self.nx + 1) * (self.ny + 1) * (self.nz + 1)
|
||||
vtk_numpoint_components = 3
|
||||
vtk_numlines = 2 * self.nx * self.ny + 2 * self.ny * self.nz + 2 * self.nx * self.nz + 3 * self.nx * self.ny * self.nz + self.nx + self.ny + self.nz
|
||||
vtk_numline_components = 2
|
||||
vtk_numline_components = 2
|
||||
|
||||
vtk_connectivity_offset = round_value((vtk_numpoints * vtk_numpoint_components * np.dtype(np.float32).itemsize) + np.dtype(np.uint32).itemsize)
|
||||
vtk_offsets_offset = round_value(vtk_connectivity_offset + (vtk_numlines * vtk_numline_components * np.dtype(np.uint32).itemsize) + np.dtype(np.uint32).itemsize)
|
||||
@@ -244,7 +244,6 @@ class GeometryView(object):
|
||||
|
||||
self.write_gprmax_info(f, G, materialsonly=True)
|
||||
|
||||
|
||||
def write_gprmax_info(self, f, G, materialsonly=False):
|
||||
"""Writes gprMax specific information relating material, source, and receiver names to numeric identifiers.
|
||||
|
||||
|
@@ -63,6 +63,7 @@ class Grid(object):
|
||||
co = round_value(float(val) / getattr(self, 'd' + coord))
|
||||
return co
|
||||
|
||||
|
||||
class FDTDGrid(Grid):
|
||||
"""Holds attributes associated with the entire grid. A convenient way for accessing regularly used parameters."""
|
||||
|
||||
@@ -155,7 +156,7 @@ def dispersion_check(G):
|
||||
if waveform.type == 'sine' or waveform.type == 'contsine':
|
||||
maxfreqs.append(4 * waveform.freq)
|
||||
|
||||
elif waveform.type =='impulse':
|
||||
elif waveform.type == 'impulse':
|
||||
pass
|
||||
|
||||
else:
|
||||
@@ -239,6 +240,7 @@ def Ix(x, y, z, Hy, Hz, G):
|
||||
Ix = G.dy * (Hy[x, y, z - 1] - Hy[x, y, z]) + G.dz * (Hz[x, y, z] - Hz[x, y - 1, z])
|
||||
return Ix
|
||||
|
||||
|
||||
def Iy(x, y, z, Hx, Hz, G):
|
||||
"""Calculates the y-component of current at a grid position.
|
||||
|
||||
@@ -256,6 +258,7 @@ def Iy(x, y, z, Hx, Hz, G):
|
||||
Iy = G.dx * (Hx[x, y, z] - Hx[x, y, z - 1]) + G.dz * (Hz[x - 1, y, z] - Hz[x, y, z])
|
||||
return Iy
|
||||
|
||||
|
||||
def Iz(x, y, z, Hx, Hy, G):
|
||||
"""Calculates the z-component of current at a grid position.
|
||||
|
||||
|
@@ -35,6 +35,7 @@ Coordinate(x=0.1, y=0.2, z=0.3)
|
||||
from collections import namedtuple
|
||||
Coordinate_tuple = namedtuple('Coordinate', ['x', 'y', 'z'])
|
||||
|
||||
|
||||
class Coordinate(Coordinate_tuple):
|
||||
"""Subclass of a namedtuple where __str__ outputs 'x y z'"""
|
||||
def __str__(self):
|
||||
@@ -62,7 +63,7 @@ def command(cmd, *parameters):
|
||||
except TypeError as e:
|
||||
# append info about cmd and parameters to the exception:
|
||||
if not e.args:
|
||||
e.args=('', )
|
||||
e.args = ('', )
|
||||
additional_info = "Creating cmd = #{} with parameters {} -> {} failed".format(cmd, parameters, filtered_list)
|
||||
e.args = e.args + (additional_info,)
|
||||
raise e
|
||||
|
@@ -46,9 +46,9 @@ def process_geometrycmds(geometry, G):
|
||||
if len(tmp) != 6:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly five parameters')
|
||||
|
||||
xs = round_value(float(tmp[1])/G.dx)
|
||||
ys = round_value(float(tmp[2])/G.dy)
|
||||
zs = round_value(float(tmp[3])/G.dz)
|
||||
xs = round_value(float(tmp[1]) / G.dx)
|
||||
ys = round_value(float(tmp[2]) / G.dy)
|
||||
zs = round_value(float(tmp[3]) / G.dz)
|
||||
geofile = tmp[4]
|
||||
matfile = tmp[5]
|
||||
|
||||
@@ -90,17 +90,16 @@ def process_geometrycmds(geometry, G):
|
||||
if G.messages:
|
||||
print('Geometry objects from file {} inserted at {:g}m, {:g}m, {:g}m, with corresponding materials file {}.'.format(geofile, xs * G.dx, ys * G.dy, zs * G.dz, matfile))
|
||||
|
||||
|
||||
elif tmp[0] == '#edge:':
|
||||
if len(tmp) != 8:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires exactly seven parameters')
|
||||
|
||||
xs = round_value(float(tmp[1])/G.dx)
|
||||
xf = round_value(float(tmp[4])/G.dx)
|
||||
ys = round_value(float(tmp[2])/G.dy)
|
||||
yf = round_value(float(tmp[5])/G.dy)
|
||||
zs = round_value(float(tmp[3])/G.dz)
|
||||
zf = round_value(float(tmp[6])/G.dz)
|
||||
xs = round_value(float(tmp[1]) / G.dx)
|
||||
xf = round_value(float(tmp[4]) / G.dx)
|
||||
ys = round_value(float(tmp[2]) / G.dy)
|
||||
yf = round_value(float(tmp[5]) / G.dy)
|
||||
zs = round_value(float(tmp[3]) / G.dz)
|
||||
zf = round_value(float(tmp[6]) / G.dz)
|
||||
|
||||
if xs < 0 or xs > G.nx:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
|
||||
@@ -150,7 +149,6 @@ def process_geometrycmds(geometry, G):
|
||||
if G.messages:
|
||||
print('Edge from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material {} created.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, tmp[7]))
|
||||
|
||||
|
||||
elif tmp[0] == '#plate:':
|
||||
if len(tmp) < 8:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least seven parameters')
|
||||
@@ -166,12 +164,12 @@ def process_geometrycmds(geometry, G):
|
||||
else:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
|
||||
|
||||
xs = round_value(float(tmp[1])/G.dx)
|
||||
xf = round_value(float(tmp[4])/G.dx)
|
||||
ys = round_value(float(tmp[2])/G.dy)
|
||||
yf = round_value(float(tmp[5])/G.dy)
|
||||
zs = round_value(float(tmp[3])/G.dz)
|
||||
zf = round_value(float(tmp[6])/G.dz)
|
||||
xs = round_value(float(tmp[1]) / G.dx)
|
||||
xf = round_value(float(tmp[4]) / G.dx)
|
||||
ys = round_value(float(tmp[2]) / G.dy)
|
||||
yf = round_value(float(tmp[5]) / G.dy)
|
||||
zs = round_value(float(tmp[3]) / G.dz)
|
||||
zf = round_value(float(tmp[6]) / G.dz)
|
||||
|
||||
if xs < 0 or xs > G.nx:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
|
||||
@@ -259,7 +257,6 @@ def process_geometrycmds(geometry, G):
|
||||
if G.messages:
|
||||
print('Plate from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material(s) {} created.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, ', '.join(materialsrequested)))
|
||||
|
||||
|
||||
elif tmp[0] == '#triangle:':
|
||||
if len(tmp) < 12:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least eleven parameters')
|
||||
@@ -286,15 +283,15 @@ def process_geometrycmds(geometry, G):
|
||||
else:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
|
||||
|
||||
x1 = round_value(float(tmp[1])/G.dx) * G.dx
|
||||
y1 = round_value(float(tmp[2])/G.dy) * G.dy
|
||||
z1 = round_value(float(tmp[3])/G.dz) * G.dz
|
||||
x2 = round_value(float(tmp[4])/G.dx) * G.dx
|
||||
y2 = round_value(float(tmp[5])/G.dy) * G.dy
|
||||
z2 = round_value(float(tmp[6])/G.dz) * G.dz
|
||||
x3 = round_value(float(tmp[7])/G.dx) * G.dx
|
||||
y3 = round_value(float(tmp[8])/G.dy) * G.dy
|
||||
z3 = round_value(float(tmp[9])/G.dz) * G.dz
|
||||
x1 = round_value(float(tmp[1]) / G.dx) * G.dx
|
||||
y1 = round_value(float(tmp[2]) / G.dy) * G.dy
|
||||
z1 = round_value(float(tmp[3]) / G.dz) * G.dz
|
||||
x2 = round_value(float(tmp[4]) / G.dx) * G.dx
|
||||
y2 = round_value(float(tmp[5]) / G.dy) * G.dy
|
||||
z2 = round_value(float(tmp[6]) / G.dz) * G.dz
|
||||
x3 = round_value(float(tmp[7]) / G.dx) * G.dx
|
||||
y3 = round_value(float(tmp[8]) / G.dy) * G.dy
|
||||
z3 = round_value(float(tmp[9]) / G.dz) * G.dz
|
||||
thickness = float(tmp[10])
|
||||
|
||||
if x1 < 0 or x2 < 0 or x3 < 0 or x1 > G.nx or x2 > G.nx or x3 > G.nx:
|
||||
@@ -379,7 +376,6 @@ def process_geometrycmds(geometry, G):
|
||||
else:
|
||||
print('Triangle with coordinates {:g}m {:g}m {:g}m, {:g}m {:g}m {:g}m, {:g}m {:g}m {:g}m of material(s) {} created.'.format(x1, y1, z1, x2, y2, z2, x3, y3, z3, ', '.join(materialsrequested)))
|
||||
|
||||
|
||||
elif tmp[0] == '#box:':
|
||||
if len(tmp) < 8:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least seven parameters')
|
||||
@@ -406,12 +402,12 @@ def process_geometrycmds(geometry, G):
|
||||
else:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
|
||||
|
||||
xs = round_value(float(tmp[1])/G.dx)
|
||||
xf = round_value(float(tmp[4])/G.dx)
|
||||
ys = round_value(float(tmp[2])/G.dy)
|
||||
yf = round_value(float(tmp[5])/G.dy)
|
||||
zs = round_value(float(tmp[3])/G.dz)
|
||||
zf = round_value(float(tmp[6])/G.dz)
|
||||
xs = round_value(float(tmp[1]) / G.dx)
|
||||
xf = round_value(float(tmp[4]) / G.dx)
|
||||
ys = round_value(float(tmp[2]) / G.dy)
|
||||
yf = round_value(float(tmp[5]) / G.dy)
|
||||
zs = round_value(float(tmp[3]) / G.dz)
|
||||
zf = round_value(float(tmp[6]) / G.dz)
|
||||
|
||||
if xs < 0 or xs > G.nx:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
|
||||
@@ -471,7 +467,6 @@ def process_geometrycmds(geometry, G):
|
||||
dielectricsmoothing = 'off'
|
||||
print('Box from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material(s) {} created, dielectric smoothing is {}.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, ', '.join(materialsrequested), dielectricsmoothing))
|
||||
|
||||
|
||||
elif tmp[0] == '#cylinder:':
|
||||
if len(tmp) < 9:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least eight parameters')
|
||||
@@ -498,12 +493,12 @@ def process_geometrycmds(geometry, G):
|
||||
else:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
|
||||
|
||||
x1 = round_value(float(tmp[1])/G.dx) * G.dx
|
||||
y1 = round_value(float(tmp[2])/G.dy) * G.dy
|
||||
z1 = round_value(float(tmp[3])/G.dz) * G.dz
|
||||
x2 = round_value(float(tmp[4])/G.dx) * G.dx
|
||||
y2 = round_value(float(tmp[5])/G.dy) * G.dy
|
||||
z2 = round_value(float(tmp[6])/G.dz) * G.dz
|
||||
x1 = round_value(float(tmp[1]) / G.dx) * G.dx
|
||||
y1 = round_value(float(tmp[2]) / G.dy) * G.dy
|
||||
z1 = round_value(float(tmp[3]) / G.dz) * G.dz
|
||||
x2 = round_value(float(tmp[4]) / G.dx) * G.dx
|
||||
y2 = round_value(float(tmp[5]) / G.dy) * G.dy
|
||||
z2 = round_value(float(tmp[6]) / G.dz) * G.dz
|
||||
r = float(tmp[7])
|
||||
|
||||
if r <= 0:
|
||||
@@ -552,7 +547,6 @@ def process_geometrycmds(geometry, G):
|
||||
dielectricsmoothing = 'off'
|
||||
print('Cylinder with face centres {:g}m, {:g}m, {:g}m and {:g}m, {:g}m, {:g}m, with radius {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(x1, y1, z1, x2, y2, z2, r, ', '.join(materialsrequested), dielectricsmoothing))
|
||||
|
||||
|
||||
elif tmp[0] == '#cylindrical_sector:':
|
||||
if len(tmp) < 10:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least nine parameters')
|
||||
@@ -586,8 +580,8 @@ def process_geometrycmds(geometry, G):
|
||||
extent2 = float(tmp[5])
|
||||
thickness = extent2 - extent1
|
||||
r = float(tmp[6])
|
||||
sectorstartangle = 2*np.pi*(float(tmp[7])/360)
|
||||
sectorangle = 2*np.pi*(float(tmp[8])/360)
|
||||
sectorstartangle = 2 * np.pi * (float(tmp[7]) / 360)
|
||||
sectorangle = 2 * np.pi * (float(tmp[8]) / 360)
|
||||
|
||||
if normal != 'x' and normal != 'y' and normal != 'z':
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the normal direction must be either x, y or z.')
|
||||
@@ -648,21 +642,21 @@ def process_geometrycmds(geometry, G):
|
||||
|
||||
# yz-plane cylindrical sector
|
||||
if normal == 'x':
|
||||
ctr1 = round_value(ctr1/G.dy) * G.dy
|
||||
ctr2 = round_value(ctr2/G.dz) * G.dz
|
||||
level = round_value(extent1/G.dx)
|
||||
ctr1 = round_value(ctr1 / G.dy) * G.dy
|
||||
ctr2 = round_value(ctr2 / G.dz) * G.dz
|
||||
level = round_value(extent1 / G.dx)
|
||||
|
||||
# xz-plane cylindrical sector
|
||||
elif normal == 'y':
|
||||
ctr1 = round_value(ctr1/G.dx) * G.dx
|
||||
ctr2 = round_value(ctr2/G.dz) * G.dz
|
||||
level = round_value(extent1/G.dy)
|
||||
ctr1 = round_value(ctr1 / G.dx) * G.dx
|
||||
ctr2 = round_value(ctr2 / G.dz) * G.dz
|
||||
level = round_value(extent1 / G.dy)
|
||||
|
||||
# xy-plane cylindrical sector
|
||||
elif normal == 'z':
|
||||
ctr1 = round_value(ctr1/G.dx) * G.dx
|
||||
ctr2 = round_value(ctr2/G.dy) * G.dy
|
||||
level = round_value(extent1/G.dz)
|
||||
ctr1 = round_value(ctr1 / G.dx) * G.dx
|
||||
ctr2 = round_value(ctr2 / G.dy) * G.dy
|
||||
level = round_value(extent1 / G.dz)
|
||||
|
||||
build_cylindrical_sector(ctr1, ctr2, level, sectorstartangle, sectorangle, r, normal, thickness, G.dx, G.dy, G.dz, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID)
|
||||
|
||||
@@ -672,10 +666,9 @@ def process_geometrycmds(geometry, G):
|
||||
dielectricsmoothing = 'on'
|
||||
else:
|
||||
dielectricsmoothing = 'off'
|
||||
print('Cylindrical sector with centre {:g}m, {:g}m, radius {:g}m, starting angle {:.1f} degrees, sector angle {:.1f} degrees, thickness {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(ctr1, ctr2, r, (sectorstartangle/(2*np.pi))*360, (sectorangle/(2*np.pi))*360, thickness, ', '.join(materialsrequested), dielectricsmoothing))
|
||||
print('Cylindrical sector with centre {:g}m, {:g}m, radius {:g}m, starting angle {:.1f} degrees, sector angle {:.1f} degrees, thickness {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(ctr1, ctr2, r, (sectorstartangle / (2 * np.pi)) * 360, (sectorangle / (2 * np.pi)) * 360, thickness, ', '.join(materialsrequested), dielectricsmoothing))
|
||||
else:
|
||||
print('Cylindrical sector with centre {:g}m, {:g}m, radius {:g}m, starting angle {:.1f} degrees, sector angle {:.1f} degrees, of material(s) {} created.'.format(ctr1, ctr2, r, (sectorstartangle/(2*np.pi))*360, (sectorangle/(2*np.pi))*360, ', '.join(materialsrequested)))
|
||||
|
||||
print('Cylindrical sector with centre {:g}m, {:g}m, radius {:g}m, starting angle {:.1f} degrees, sector angle {:.1f} degrees, of material(s) {} created.'.format(ctr1, ctr2, r, (sectorstartangle / (2 * np.pi)) * 360, (sectorangle / (2 * np.pi)) * 360, ', '.join(materialsrequested)))
|
||||
|
||||
elif tmp[0] == '#sphere:':
|
||||
if len(tmp) < 6:
|
||||
@@ -704,9 +697,9 @@ def process_geometrycmds(geometry, G):
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
|
||||
|
||||
# Centre of sphere
|
||||
xc = round_value(float(tmp[1])/G.dx)
|
||||
yc = round_value(float(tmp[2])/G.dy)
|
||||
zc = round_value(float(tmp[3])/G.dz)
|
||||
xc = round_value(float(tmp[1]) / G.dx)
|
||||
yc = round_value(float(tmp[2]) / G.dy)
|
||||
zc = round_value(float(tmp[3]) / G.dz)
|
||||
r = float(tmp[4])
|
||||
|
||||
# Look up requested materials in existing list of material instances
|
||||
@@ -752,7 +745,6 @@ def process_geometrycmds(geometry, G):
|
||||
dielectricsmoothing = 'off'
|
||||
print('Sphere with centre {:g}m, {:g}m, {:g}m, radius {:g}m, of material(s) {} created, dielectric smoothing is {}.'.format(xc * G.dx, yc * G.dy, zc * G.dz, r, ', '.join(materialsrequested), dielectricsmoothing))
|
||||
|
||||
|
||||
elif tmp[0] == '#fractal_box:':
|
||||
if len(tmp) < 14:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least thirteen parameters')
|
||||
@@ -763,12 +755,12 @@ def process_geometrycmds(geometry, G):
|
||||
else:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' too many parameters have been given')
|
||||
|
||||
xs = round_value(float(tmp[1])/G.dx)
|
||||
xf = round_value(float(tmp[4])/G.dx)
|
||||
ys = round_value(float(tmp[2])/G.dy)
|
||||
yf = round_value(float(tmp[5])/G.dy)
|
||||
zs = round_value(float(tmp[3])/G.dz)
|
||||
zf = round_value(float(tmp[6])/G.dz)
|
||||
xs = round_value(float(tmp[1]) / G.dx)
|
||||
xf = round_value(float(tmp[4]) / G.dx)
|
||||
ys = round_value(float(tmp[2]) / G.dy)
|
||||
yf = round_value(float(tmp[5]) / G.dy)
|
||||
zs = round_value(float(tmp[3]) / G.dz)
|
||||
zf = round_value(float(tmp[6]) / G.dz)
|
||||
|
||||
if xs < 0 or xs > G.nx:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
|
||||
@@ -836,12 +828,12 @@ def process_geometrycmds(geometry, G):
|
||||
|
||||
# Only process rough surfaces for this fractal volume
|
||||
if tmp[12] == volume.ID:
|
||||
xs = round_value(float(tmp[1])/G.dx)
|
||||
xf = round_value(float(tmp[4])/G.dx)
|
||||
ys = round_value(float(tmp[2])/G.dy)
|
||||
yf = round_value(float(tmp[5])/G.dy)
|
||||
zs = round_value(float(tmp[3])/G.dz)
|
||||
zf = round_value(float(tmp[6])/G.dz)
|
||||
xs = round_value(float(tmp[1]) / G.dx)
|
||||
xf = round_value(float(tmp[4]) / G.dx)
|
||||
ys = round_value(float(tmp[2]) / G.dy)
|
||||
yf = round_value(float(tmp[5]) / G.dy)
|
||||
zs = round_value(float(tmp[3]) / G.dz)
|
||||
zf = round_value(float(tmp[6]) / G.dz)
|
||||
|
||||
if xs < 0 or xs > G.nx:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
|
||||
@@ -943,12 +935,12 @@ def process_geometrycmds(geometry, G):
|
||||
|
||||
# Only process surfaces for this fractal volume
|
||||
if tmp[8] == volume.ID:
|
||||
xs = round_value(float(tmp[1])/G.dx)
|
||||
xf = round_value(float(tmp[4])/G.dx)
|
||||
ys = round_value(float(tmp[2])/G.dy)
|
||||
yf = round_value(float(tmp[5])/G.dy)
|
||||
zs = round_value(float(tmp[3])/G.dz)
|
||||
zf = round_value(float(tmp[6])/G.dz)
|
||||
xs = round_value(float(tmp[1]) / G.dx)
|
||||
xf = round_value(float(tmp[4]) / G.dx)
|
||||
ys = round_value(float(tmp[2]) / G.dy)
|
||||
yf = round_value(float(tmp[5]) / G.dy)
|
||||
zs = round_value(float(tmp[3]) / G.dz)
|
||||
zf = round_value(float(tmp[6]) / G.dz)
|
||||
depth = float(tmp[7])
|
||||
|
||||
if xs < 0 or xs > G.nx:
|
||||
@@ -1044,7 +1036,6 @@ def process_geometrycmds(geometry, G):
|
||||
if G.messages:
|
||||
print('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 * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, filldepth, surface.operatingonID))
|
||||
|
||||
|
||||
if tmp[0] == '#add_grass:':
|
||||
if len(tmp) < 12:
|
||||
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least eleven parameters')
|
||||
@@ -1057,12 +1048,12 @@ def process_geometrycmds(geometry, G):
|
||||
|
||||
# Only process grass for this fractal volume
|
||||
if tmp[11] == volume.ID:
|
||||
xs = round_value(float(tmp[1])/G.dx)
|
||||
xf = round_value(float(tmp[4])/G.dx)
|
||||
ys = round_value(float(tmp[2])/G.dy)
|
||||
yf = round_value(float(tmp[5])/G.dy)
|
||||
zs = round_value(float(tmp[3])/G.dz)
|
||||
zf = round_value(float(tmp[6])/G.dz)
|
||||
xs = round_value(float(tmp[1]) / G.dx)
|
||||
xf = round_value(float(tmp[4]) / G.dx)
|
||||
ys = round_value(float(tmp[2]) / G.dy)
|
||||
yf = round_value(float(tmp[5]) / G.dy)
|
||||
zs = round_value(float(tmp[3]) / G.dz)
|
||||
zf = round_value(float(tmp[6]) / G.dz)
|
||||
numblades = int(tmp[10])
|
||||
|
||||
if xs < 0 or xs > G.nx:
|
||||
@@ -1194,7 +1185,6 @@ def process_geometrycmds(geometry, G):
|
||||
if G.messages:
|
||||
print('{} 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(numblades, xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, surface.dimension, surface.seed, float(tmp[8]), float(tmp[9]), surface.operatingonID))
|
||||
|
||||
|
||||
# Process any modifications to the original fractal box then generate it
|
||||
if volume.fractalsurfaces:
|
||||
volume.originalxs = volume.xs
|
||||
@@ -1439,7 +1429,6 @@ def process_geometrycmds(geometry, G):
|
||||
depth += 1
|
||||
k -= 1
|
||||
|
||||
|
||||
# Build voxels from any true values of the 3D mask array
|
||||
waternumID = next((x.numID for x in G.materials if x.ID == 'water'), 0)
|
||||
grassnumID = next((x.numID for x in G.materials if x.ID == 'grass'), 0)
|
||||
|
@@ -71,7 +71,6 @@ def process_multicmds(multicmds, G):
|
||||
|
||||
G.waveforms.append(w)
|
||||
|
||||
|
||||
# Voltage source
|
||||
cmdname = '#voltage_source'
|
||||
if multicmds[cmdname] != 'None':
|
||||
@@ -98,7 +97,7 @@ def process_multicmds(multicmds, G):
|
||||
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' there is no waveform with the identifier {}'.format(tmp[5]))
|
||||
|
||||
v = VoltageSource()
|
||||
v.polarisation= tmp[0]
|
||||
v.polarisation = tmp[0]
|
||||
v.xcoord = xcoord
|
||||
v.ycoord = ycoord
|
||||
v.zcoord = zcoord
|
||||
@@ -132,7 +131,6 @@ def process_multicmds(multicmds, G):
|
||||
|
||||
G.voltagesources.append(v)
|
||||
|
||||
|
||||
# Hertzian dipole
|
||||
cmdname = '#hertzian_dipole'
|
||||
if multicmds[cmdname] != 'None':
|
||||
@@ -192,7 +190,6 @@ def process_multicmds(multicmds, G):
|
||||
|
||||
G.hertziandipoles.append(h)
|
||||
|
||||
|
||||
# Magnetic dipole
|
||||
cmdname = '#magnetic_dipole'
|
||||
if multicmds[cmdname] != 'None':
|
||||
@@ -325,9 +322,9 @@ def process_multicmds(multicmds, G):
|
||||
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' has an incorrect number of parameters')
|
||||
|
||||
# Check position parameters
|
||||
xcoord = round_value(float(tmp[0])/G.dx)
|
||||
ycoord = round_value(float(tmp[1])/G.dy)
|
||||
zcoord = round_value(float(tmp[2])/G.dz)
|
||||
xcoord = round_value(float(tmp[0]) / G.dx)
|
||||
ycoord = round_value(float(tmp[1]) / G.dy)
|
||||
zcoord = round_value(float(tmp[2]) / G.dz)
|
||||
check_coordinates(xcoord, ycoord, zcoord)
|
||||
if xcoord < G.pmlthickness[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]:
|
||||
print("WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.\n')
|
||||
@@ -359,7 +356,6 @@ def process_multicmds(multicmds, G):
|
||||
|
||||
G.rxs.append(r)
|
||||
|
||||
|
||||
# Receiver array
|
||||
cmdname = '#rx_array'
|
||||
if multicmds[cmdname] != 'None':
|
||||
@@ -425,7 +421,6 @@ def process_multicmds(multicmds, G):
|
||||
print(' Receiver at {:g}m, {:g}m, {:g}m with output(s) {} created.'.format(r.xcoord * G.dx, r.ycoord * G.dy, r.zcoord * G.dz, ', '.join(r.outputs)))
|
||||
G.rxs.append(r)
|
||||
|
||||
|
||||
# Snapshot
|
||||
cmdname = '#snapshot'
|
||||
if multicmds[cmdname] != 'None':
|
||||
@@ -511,7 +506,6 @@ def process_multicmds(multicmds, G):
|
||||
# Append the new material object to the materials list
|
||||
G.materials.append(m)
|
||||
|
||||
|
||||
cmdname = '#add_dispersion_debye'
|
||||
if multicmds[cmdname] != 'None':
|
||||
for cmdinstance in multicmds[cmdname]:
|
||||
@@ -583,7 +577,6 @@ def process_multicmds(multicmds, G):
|
||||
if G.messages:
|
||||
print('Lorentz disperion added to {} with delta_epsr={}, 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)))
|
||||
|
||||
|
||||
cmdname = '#add_dispersion_drude'
|
||||
if multicmds[cmdname] != 'None':
|
||||
for cmdinstance in multicmds[cmdname]:
|
||||
@@ -619,7 +612,6 @@ def process_multicmds(multicmds, G):
|
||||
if G.messages:
|
||||
print('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)))
|
||||
|
||||
|
||||
cmdname = '#soil_peplinski'
|
||||
if multicmds[cmdname] != 'None':
|
||||
for cmdinstance in multicmds[cmdname]:
|
||||
@@ -650,7 +642,6 @@ def process_multicmds(multicmds, G):
|
||||
# Append the new material object to the materials list
|
||||
G.mixingmodels.append(s)
|
||||
|
||||
|
||||
# Geometry views (creates VTK-based geometry files)
|
||||
cmdname = '#geometry_view'
|
||||
if multicmds[cmdname] != 'None':
|
||||
@@ -695,7 +686,6 @@ def process_multicmds(multicmds, G):
|
||||
# Append the new GeometryView object to the geometry views list
|
||||
G.geometryviews.append(g)
|
||||
|
||||
|
||||
# Complex frequency shifted (CFS) PML parameter
|
||||
cmdname = '#pml_cfs'
|
||||
if multicmds[cmdname] != 'None':
|
||||
|
@@ -50,7 +50,6 @@ def process_singlecmds(singlecmds, G):
|
||||
else:
|
||||
raise CmdInputError(cmd + ' requires input values of either y or n')
|
||||
|
||||
|
||||
# Title
|
||||
cmd = '#title'
|
||||
if singlecmds[cmd] != 'None':
|
||||
@@ -58,7 +57,6 @@ def process_singlecmds(singlecmds, G):
|
||||
if G.messages:
|
||||
print('Model title: {}'.format(G.title))
|
||||
|
||||
|
||||
# Number of processors to run on (OpenMP)
|
||||
cmd = '#num_threads'
|
||||
os.environ['OMP_WAIT_POLICY'] = 'active'
|
||||
@@ -85,7 +83,6 @@ def process_singlecmds(singlecmds, G):
|
||||
if G.nthreads > psutil.cpu_count(logical=False):
|
||||
print('\nWARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, psutil.cpu_count(logical=False)))
|
||||
|
||||
|
||||
# Spatial discretisation
|
||||
cmd = '#dx_dy_dz'
|
||||
tmp = [float(x) for x in singlecmds[cmd].split()]
|
||||
@@ -103,15 +100,14 @@ def process_singlecmds(singlecmds, G):
|
||||
if G.messages:
|
||||
print('Spatial discretisation: {:g} x {:g} x {:g}m'.format(G.dx, G.dy, G.dz))
|
||||
|
||||
|
||||
# Domain
|
||||
cmd = '#domain'
|
||||
tmp = [float(x) for x in singlecmds[cmd].split()]
|
||||
if len(tmp) != 3:
|
||||
raise CmdInputError(cmd + ' requires exactly three parameters')
|
||||
G.nx = round_value(tmp[0]/G.dx)
|
||||
G.ny = round_value(tmp[1]/G.dy)
|
||||
G.nz = round_value(tmp[2]/G.dz)
|
||||
G.nx = round_value(tmp[0] / G.dx)
|
||||
G.ny = round_value(tmp[1] / G.dy)
|
||||
G.nz = round_value(tmp[2] / G.dz)
|
||||
if G.nx == 0 or G.ny == 0 or G.nz == 0:
|
||||
raise CmdInputError(cmd + ' requires at least one cell in every dimension')
|
||||
if G.messages:
|
||||
@@ -120,7 +116,6 @@ def process_singlecmds(singlecmds, G):
|
||||
mem = (((G.nx + 1) * (G.ny + 1) * (G.nz + 1) * 13 * np.dtype(floattype).itemsize + (G.nx + 1) * (G.ny + 1) * (G.nz + 1) * 18) * 1.1) + 30e6
|
||||
print('Memory (RAM) usage: ~{} required, {} available'.format(human_size(mem), human_size(psutil.virtual_memory().total)))
|
||||
|
||||
|
||||
# Time step CFL limit (use either 2D or 3D) and default PML thickness
|
||||
if G.nx == 1:
|
||||
G.dt = 1 / (c * np.sqrt((1 / G.dy) * (1 / G.dy) + (1 / G.dz) * (1 / G.dz)))
|
||||
@@ -145,7 +140,6 @@ def process_singlecmds(singlecmds, G):
|
||||
if G.messages:
|
||||
print('Time step (at {} CFL limit): {:g} secs'.format(G.dtlimit, G.dt))
|
||||
|
||||
|
||||
# Time step stability factor
|
||||
cmd = '#time_step_stability_factor'
|
||||
if singlecmds[cmd] != 'None':
|
||||
@@ -158,7 +152,6 @@ def process_singlecmds(singlecmds, G):
|
||||
if G.messages:
|
||||
print('Time step (modified): {:g} secs'.format(G.dt))
|
||||
|
||||
|
||||
# Time window
|
||||
cmd = '#time_window'
|
||||
tmp = singlecmds[cmd].split()
|
||||
@@ -182,7 +175,6 @@ def process_singlecmds(singlecmds, G):
|
||||
if G.messages:
|
||||
print('Time window: {:g} secs ({} iterations)'.format(G.timewindow, G.iterations))
|
||||
|
||||
|
||||
# PML
|
||||
cmd = '#pml_cells'
|
||||
if singlecmds[cmd] != 'None':
|
||||
@@ -193,36 +185,33 @@ def process_singlecmds(singlecmds, G):
|
||||
G.pmlthickness = (int(tmp[0]), int(tmp[0]), int(tmp[0]), int(tmp[0]), int(tmp[0]), int(tmp[0]))
|
||||
else:
|
||||
G.pmlthickness = (int(tmp[0]), int(tmp[1]), int(tmp[2]), int(tmp[3]), int(tmp[4]), int(tmp[5]))
|
||||
if 2*G.pmlthickness[0] >= G.nx or 2*G.pmlthickness[1] >= G.ny or 2*G.pmlthickness[2] >= G.nz or 2*G.pmlthickness[3] >= G.nx or 2*G.pmlthickness[4] >= G.ny or 2*G.pmlthickness[5] >= G.nz:
|
||||
if 2 * G.pmlthickness[0] >= G.nx or 2 * G.pmlthickness[1] >= G.ny or 2 * G.pmlthickness[2] >= G.nz or 2 * G.pmlthickness[3] >= G.nx or 2 * G.pmlthickness[4] >= G.ny or 2 * G.pmlthickness[5] >= G.nz:
|
||||
raise CmdInputError(cmd + ' has too many cells for the domain size')
|
||||
|
||||
|
||||
# src_steps
|
||||
cmd = '#src_steps'
|
||||
if singlecmds[cmd] != 'None':
|
||||
tmp = singlecmds[cmd].split()
|
||||
if len(tmp) != 3:
|
||||
raise CmdInputError(cmd + ' requires exactly three parameters')
|
||||
G.srcstepx = round_value(float(tmp[0])/G.dx)
|
||||
G.srcstepy = round_value(float(tmp[1])/G.dy)
|
||||
G.srcstepz = round_value(float(tmp[2])/G.dz)
|
||||
G.srcstepx = round_value(float(tmp[0]) / G.dx)
|
||||
G.srcstepy = round_value(float(tmp[1]) / G.dy)
|
||||
G.srcstepz = round_value(float(tmp[2]) / G.dz)
|
||||
if G.messages:
|
||||
print('Simple sources will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.srcstepx * G.dx, G.srcstepy * G.dy, G.srcstepz * G.dz))
|
||||
|
||||
|
||||
# rx_steps
|
||||
cmd = '#rx_steps'
|
||||
if singlecmds[cmd] != 'None':
|
||||
tmp = singlecmds[cmd].split()
|
||||
if len(tmp) != 3:
|
||||
raise CmdInputError(cmd + ' requires exactly three parameters')
|
||||
G.rxstepx = round_value(float(tmp[0])/G.dx)
|
||||
G.rxstepy = round_value(float(tmp[1])/G.dy)
|
||||
G.rxstepz = round_value(float(tmp[2])/G.dz)
|
||||
G.rxstepx = round_value(float(tmp[0]) / G.dx)
|
||||
G.rxstepy = round_value(float(tmp[1]) / G.dy)
|
||||
G.rxstepz = round_value(float(tmp[2]) / G.dz)
|
||||
if G.messages:
|
||||
print('All receivers will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.rxstepx * G.dx, G.rxstepy * G.dy, G.rxstepz * G.dz))
|
||||
|
||||
|
||||
# Excitation file for user-defined source waveforms
|
||||
cmd = '#excitation_file'
|
||||
if singlecmds[cmd] != 'None':
|
||||
@@ -251,7 +240,7 @@ def process_singlecmds(singlecmds, G):
|
||||
if len(waveformvalues.shape) == 1:
|
||||
w.uservalues = waveformvalues[:]
|
||||
else:
|
||||
w.uservalues = waveformvalues[:,waveform]
|
||||
w.uservalues = waveformvalues[:, waveform]
|
||||
|
||||
if G.messages:
|
||||
print('User waveform {} created.'.format(w.ID))
|
||||
|
@@ -74,8 +74,8 @@ class Material(object):
|
||||
G (class): Grid class instance - holds essential parameters describing the model.
|
||||
"""
|
||||
|
||||
HA = (m0*self.mr / G.dt) + 0.5*self.sm
|
||||
HB = (m0*self.mr / G.dt) - 0.5*self.sm
|
||||
HA = (m0 * self.mr / G.dt) + 0.5 * self.sm
|
||||
HB = (m0 * self.mr / G.dt) - 0.5 * self.sm
|
||||
self.DA = HB / HA
|
||||
self.DBx = (1 / G.dx) * 1 / HA
|
||||
self.DBy = (1 / G.dy) * 1 / HA
|
||||
@@ -122,12 +122,12 @@ class Material(object):
|
||||
self.zt[x] = (self.w[x] / self.q[x]) * (1 - self.eqt[x]) / G.dt
|
||||
self.zt2[x] = (self.w[x] / self.q[x]) * (1 - self.eqt2[x])
|
||||
|
||||
EA = (e0*self.er / G.dt) + 0.5*self.se - (e0 / G.dt) * np.sum(self.zt2.real)
|
||||
EB = (e0*self.er / G.dt) - 0.5*self.se - (e0 / G.dt) * np.sum(self.zt2.real)
|
||||
EA = (e0 * self.er / G.dt) + 0.5 * self.se - (e0 / G.dt) * np.sum(self.zt2.real)
|
||||
EB = (e0 * self.er / G.dt) - 0.5 * self.se - (e0 / G.dt) * np.sum(self.zt2.real)
|
||||
|
||||
else:
|
||||
EA = (e0*self.er / G.dt) + 0.5*self.se
|
||||
EB = (e0*self.er / G.dt) - 0.5*self.se
|
||||
EA = (e0 * self.er / G.dt) + 0.5 * self.se
|
||||
EB = (e0 * self.er / G.dt) - 0.5 * self.se
|
||||
|
||||
if self.ID == 'pec' or self.se == float('inf'):
|
||||
self.CA = 0
|
||||
@@ -198,15 +198,15 @@ class PeplinskiSoil(object):
|
||||
muiter = np.nditer(mumaterials, flags=['c_index'])
|
||||
while not muiter.finished:
|
||||
# Real part for frequencies in the range 1.4GHz to 18GHz
|
||||
er1 = (1 + (self.rb/self.rs) * ((es**a) - 1) + (muiter[0]**b1 * erealw**a) - muiter[0]) ** (1/a)
|
||||
er1 = (1 + (self.rb / self.rs) * ((es**a) - 1) + (muiter[0]**b1 * erealw**a) - muiter[0]) ** (1 / a)
|
||||
# Real part for frequencies in the range 0.3GHz to 1.3GHz
|
||||
er2 = 1.15 * er1 - 0.68
|
||||
|
||||
# Imaginary part for frequencies in the range 0.3GHz to 1.3GHz
|
||||
eri = er2 - (muiter[0]**(b2/a) * Material.waterdeltaer)
|
||||
eri = er2 - (muiter[0]**(b2 / a) * Material.waterdeltaer)
|
||||
|
||||
# Effective conductivity
|
||||
sig = muiter[0]**(b2/a) * ((sigf1 * (self.rs - self.rb)) / (self.rs * muiter[0]))
|
||||
sig = muiter[0]**(b2 / a) * ((sigf1 * (self.rs - self.rb)) / (self.rs * muiter[0]))
|
||||
|
||||
# Check to see if the material already exists before creating a new one
|
||||
requiredID = '|{:.4f}|'.format(float(muiter[0]))
|
||||
|
@@ -68,7 +68,7 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace):
|
||||
|
||||
# Select OA
|
||||
OA, N, cols, k, s, t = construct_OA(optparams)
|
||||
print('\n{}\nTaguchi optimisation...\n'.format(68*'*'))
|
||||
print('\n{}\nTaguchi optimisation...\n'.format(68 * '*'))
|
||||
print('\tOrthogonal array: {:g} experiments per iteration, {:g} parameters ({:g} will be used), {:g} levels, and strength {:g}'.format(N, cols, k, s, t))
|
||||
tmp = [(k, v) for k, v in optparams.items()]
|
||||
print('\tParameters to optimise with ranges: {}'.format(str(tmp).strip('[]')))
|
||||
@@ -131,7 +131,7 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace):
|
||||
# Rename confirmation experiment output file so that it is retained for each iteraction
|
||||
os.rename(outputfile, os.path.splitext(outputfile)[0] + '_final' + str(iteration + 1) + '.out')
|
||||
|
||||
print('\nTaguchi optimisation, iteration {} completed. History of optimal parameter values {} and of fitness values {}'.format(iteration + 1, dict(optparamshist), fitnessvalueshist, 68*'*'))
|
||||
print('\nTaguchi optimisation, iteration {} completed. History of optimal parameter values {} and of fitness values {}'.format(iteration + 1, dict(optparamshist), fitnessvalueshist, 68 * '*'))
|
||||
iteration += 1
|
||||
|
||||
# Stop optimisation if stopping criterion has been reached
|
||||
@@ -154,7 +154,7 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace):
|
||||
pickle.dump(fitnessvalueshist, f)
|
||||
pickle.dump(optparamsinit, f)
|
||||
|
||||
print('\n{}\nTaguchi optimisation completed after {} iteration(s).\nHistory of optimal parameter values {} and of fitness values {}\n{}\n'.format(68*'*', iteration, dict(optparamshist), fitnessvalueshist, 68*'*'))
|
||||
print('\n{}\nTaguchi optimisation completed after {} iteration(s).\nHistory of optimal parameter values {} and of fitness values {}\n{}\n'.format(68 * '*', iteration, dict(optparamshist), fitnessvalueshist, 68 * '*'))
|
||||
|
||||
|
||||
def taguchi_code_blocks(inputfile, taguchinamespace):
|
||||
@@ -313,7 +313,7 @@ def calculate_ranges_experiments(optparams, optparamsinit, levels, levelsopt, le
|
||||
|
||||
# Gaussian reduction function used for calculating levels
|
||||
T = 18 # Usually values between 15 - 20
|
||||
RR = np.exp(-(i/T)**2)
|
||||
RR = np.exp(-(i / T)**2)
|
||||
|
||||
# Calculate levels for each parameter
|
||||
for p in range(k):
|
||||
@@ -341,7 +341,6 @@ def calculate_ranges_experiments(optparams, optparamsinit, levels, levelsopt, le
|
||||
levels[0, p] = levels[1, p] - levelsdiff[p]
|
||||
levels[2, p] = levels[1, p] + levelsdiff[p]
|
||||
|
||||
|
||||
# Update dictionary of parameters to optimise with lists of new values; clear dictionary first
|
||||
optparams = OrderedDict((key, list()) for key in optparams)
|
||||
p = 0
|
||||
|
@@ -98,7 +98,7 @@ class CFS(object):
|
||||
Hvalues (float): numpy array holding scaling profile values for magnetic PML update.
|
||||
"""
|
||||
|
||||
tmp = (np.linspace(0, (len(Evalues) - 1) + 0.5, num=2*len(Evalues)) / (len(Evalues) - 1)) ** order
|
||||
tmp = (np.linspace(0, (len(Evalues) - 1) + 0.5, num=2 * len(Evalues)) / (len(Evalues) - 1)) ** order
|
||||
Evalues = tmp[0:-1:2]
|
||||
Hvalues = tmp[1::2]
|
||||
return Evalues, Hvalues
|
||||
@@ -115,7 +115,7 @@ class CFS(object):
|
||||
Hvalues (float): numpy array holding profile value for magnetic PML update.
|
||||
"""
|
||||
|
||||
Evalues= np.zeros(thickness + 1, dtype=floattype)
|
||||
Evalues = np.zeros(thickness + 1, dtype=floattype)
|
||||
Hvalues = np.zeros(thickness + 1, dtype=floattype)
|
||||
|
||||
if parameter.scalingprofile == 'constant':
|
||||
@@ -218,18 +218,18 @@ class PML(object):
|
||||
Esigma, Hsigma = cfs.calculate_values(self.thickness, cfs.sigma)
|
||||
|
||||
# Electric PML update coefficients
|
||||
tmp = (2*e0*Ekappa) + G.dt * (Ealpha * Ekappa + Esigma)
|
||||
self.ERA[x, :] = (2*e0 + G.dt*Ealpha) / tmp
|
||||
self.ERB[x, :] = (2*e0*Ekappa) / tmp
|
||||
self.ERE[x, :] = ((2*e0*Ekappa) - G.dt * (Ealpha * Ekappa + Esigma)) / tmp
|
||||
self.ERF[x, :] = (2*Esigma*G.dt) / (Ekappa * tmp)
|
||||
tmp = (2 * e0 * Ekappa) + G.dt * (Ealpha * Ekappa + Esigma)
|
||||
self.ERA[x, :] = (2 * e0 + G.dt * Ealpha) / tmp
|
||||
self.ERB[x, :] = (2 * e0 * Ekappa) / tmp
|
||||
self.ERE[x, :] = ((2 * e0 * Ekappa) - G.dt * (Ealpha * Ekappa + Esigma)) / tmp
|
||||
self.ERF[x, :] = (2 * Esigma * G.dt) / (Ekappa * tmp)
|
||||
|
||||
# Magnetic PML update coefficients
|
||||
tmp = (2*e0*Hkappa) + G.dt * (Halpha * Hkappa + Hsigma)
|
||||
self.HRA[x, :] = (2*e0 + G.dt*Halpha) / tmp
|
||||
self.HRB[x, :] = (2*e0*Hkappa) / tmp
|
||||
self.HRE[x, :] = ((2*e0*Hkappa) - G.dt * (Halpha * Hkappa + Hsigma)) / tmp
|
||||
self.HRF[x, :] = (2*Hsigma*G.dt) / (Hkappa * tmp)
|
||||
tmp = (2 * e0 * Hkappa) + G.dt * (Halpha * Hkappa + Hsigma)
|
||||
self.HRA[x, :] = (2 * e0 + G.dt * Halpha) / tmp
|
||||
self.HRB[x, :] = (2 * e0 * Hkappa) / tmp
|
||||
self.HRE[x, :] = ((2 * e0 * Hkappa) - G.dt * (Halpha * Hkappa + Hsigma)) / tmp
|
||||
self.HRF[x, :] = (2 * Hsigma * G.dt) / (Hkappa * tmp)
|
||||
|
||||
|
||||
def build_pmls(G):
|
||||
@@ -244,7 +244,7 @@ def build_pmls(G):
|
||||
if pmlthickness > 0:
|
||||
sumer = 0 # Sum of relative permittivities in PML slab
|
||||
summr = 0 # Sum of relative permeabilities in PML slab
|
||||
pmldirection=PML.directions[index]
|
||||
pmldirection = PML.directions[index]
|
||||
|
||||
if pmldirection[0] == 'x':
|
||||
if pmldirection == 'xminus':
|
||||
|
@@ -185,13 +185,13 @@ class MagneticDipole(Source):
|
||||
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
|
||||
|
||||
if self.polarisation is 'x':
|
||||
Hx[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz))
|
||||
Hx[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz))
|
||||
|
||||
elif self.polarisation is 'y':
|
||||
Hy[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz))
|
||||
Hy[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz))
|
||||
|
||||
elif self.polarisation is 'z':
|
||||
Hz[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz))
|
||||
Hz[i, j, k] -= waveform.amp * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz))
|
||||
|
||||
|
||||
class TransmissionLine(Source):
|
||||
|
@@ -57,9 +57,9 @@ along with gprMax. If not, see <http://www.gnu.org/licenses/>."""
|
||||
url = 'www.gprmax.com'
|
||||
|
||||
print('\nElectromagnetic modelling software based on the Finite-Difference \nTime-Domain (FDTD) method')
|
||||
print('\n{} {} {}'.format('*'*round((width - len(url))/2), url, '*'*round((width - len(url))/2)))
|
||||
print('\n{} {} {}'.format('*' * round((width - len(url)) / 2), url, '*' * round((width - len(url)) / 2)))
|
||||
print('{}'.format(gprMaxlogo))
|
||||
print('{} v{} {}'.format('*'*round((width - len(version))/2), (version), '*'*round((width - len(version))/2)))
|
||||
print('{} v{} {}'.format('*' * round((width - len(version)) / 2), (version), '*' * round((width - len(version)) / 2)))
|
||||
print(licenseinfo)
|
||||
|
||||
|
||||
@@ -91,12 +91,12 @@ def round_value(value, decimalplaces=0):
|
||||
|
||||
# Rounds to nearest integer (half values are rounded downwards)
|
||||
if decimalplaces == 0:
|
||||
rounded = int(d.Decimal(value).quantize(d.Decimal('1'),rounding=d.ROUND_HALF_DOWN))
|
||||
rounded = int(d.Decimal(value).quantize(d.Decimal('1'), rounding=d.ROUND_HALF_DOWN))
|
||||
|
||||
# Rounds down to nearest float represented by number of decimal places
|
||||
else:
|
||||
precision = '1.{places}'.format(places='0' * decimalplaces)
|
||||
rounded = float(d.Decimal(value).quantize(d.Decimal(precision),rounding=d.ROUND_FLOOR))
|
||||
rounded = float(d.Decimal(value).quantize(d.Decimal(precision), rounding=d.ROUND_FLOOR))
|
||||
|
||||
return rounded
|
||||
|
||||
|
@@ -43,7 +43,6 @@ def write_hdf5(outputfile, Ex, Ey, Ez, Hx, Hy, Hz, G):
|
||||
f.attrs['srcsteps'] = (G.srcstepx, G.srcstepy, G.srcstepz)
|
||||
f.attrs['rxsteps'] = (G.rxstepx, G.rxstepy, G.rxstepz)
|
||||
|
||||
|
||||
# Create group for sources (except transmission lines); add type and positional data attributes
|
||||
srclist = G.voltagesources + G.hertziandipoles + G.magneticdipoles
|
||||
for srcindex, src in enumerate(srclist):
|
||||
@@ -65,7 +64,6 @@ def write_hdf5(outputfile, Ex, Ey, Ez, Hx, Hy, Hz, G):
|
||||
f['/tls/tl' + str(tlindex + 1) + '/Vtotal'] = tl.Vtotal
|
||||
f['/tls/tl' + str(tlindex + 1) + '/Itotal'] = tl.Itotal
|
||||
|
||||
|
||||
# Create group, add positional data and write field component arrays for receivers
|
||||
for rxindex, rx in enumerate(G.rxs):
|
||||
grp = f.create_group('/rxs/rx' + str(rxindex + 1))
|
||||
|
@@ -6,6 +6,7 @@ import numpy as np
|
||||
|
||||
from gprMax.grid import Grid
|
||||
|
||||
|
||||
class ListCounter():
|
||||
|
||||
def __init__(self, npArray):
|
||||
@@ -16,10 +17,10 @@ class ListCounter():
|
||||
self.array[self.count] = item
|
||||
self.count += 1
|
||||
|
||||
|
||||
class EdgeLabels:
|
||||
|
||||
def __init__(self, grid):
|
||||
|
||||
"""
|
||||
Class to define some connectivity of for an n x l x m
|
||||
grid
|
||||
@@ -37,6 +38,7 @@ class EdgeLabels:
|
||||
edge = np.array([in_label, out_label])
|
||||
self.edge_counter.add(edge)
|
||||
|
||||
|
||||
class EdgeMaterials:
|
||||
|
||||
def __init__(self, fdtd_grid):
|
||||
@@ -85,7 +87,6 @@ class SolidLabels():
|
||||
self.label_counter = ListCounter(self.solid_labels)
|
||||
|
||||
def hexCellPicker(self, grid, i, j, k):
|
||||
|
||||
"""
|
||||
This is the ordering of nodes in the hexahedron cell.
|
||||
|
||||
@@ -124,6 +125,7 @@ class SolidLabels():
|
||||
solid_labels = self.hexCellPicker(self.label_grid.grid, i, j, k)
|
||||
self.label_counter.add(solid_labels)
|
||||
|
||||
|
||||
class SolidManager(Grid):
|
||||
|
||||
def __init__(self, label_grid, fdtd_grid):
|
||||
|
在新工单中引用
屏蔽一个用户