publicImage/gprMax
1
0
派生 0
镜像自地址 https://gitee.com/sunhf/gprMax.git 已同步 2025-08-08 07:24:19 +08:00
代码 工单 软件包 项目 版本发布 百科 动态

Merge branch 'refactor-create-internal-objects' into mpi

浏览代码
这个提交包含在：
nmannall
2024-05-21 14:45:45 +01:00
父节点 c23a2c2a40 b8b17a83bc
当前提交 818418c952
共有 32 个文件被更改，包括 1377 次插入 和 756 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

23
gprMax/cmds_geometry/add_grass.py
查看文件

@@ -103,7 +103,9 @@ class AddGrass(UserObjectGeometry):
xf, yf, zf = p2
if frac_dim < 0:
logger.exception(f"{self.__str__()} requires a positive value for the fractal dimension")
logger.exception(
f"{self.__str__()} requires a positive value for the fractal dimension"
)
raise ValueError
if limits[0] < 0 or limits[1] < 0:
logger.exception(
@@ -117,7 +119,9 @@ class AddGrass(UserObjectGeometry):
logger.exception(f"{self.__str__()} dimensions are not specified correctly")
raise ValueError
if xs not in [volume.xs, volume.xf]:
logger.exception(f"{self.__str__()} must specify external surfaces on a fractal box")
logger.exception(
f"{self.__str__()} must specify external surfaces on a fractal box"
)
raise ValueError
fractalrange = (round_value(limits[0] / grid.dx), round_value(limits[1] / grid.dx))
# xminus surface
@@ -142,7 +146,9 @@ class AddGrass(UserObjectGeometry):
logger.exception(f"{self.__str__()} dimensions are not specified correctly")
raise ValueError
if ys not in [volume.ys, volume.yf]:
logger.exception(f"{self.__str__()} must specify external surfaces on a fractal box")
logger.exception(
f"{self.__str__()} must specify external surfaces on a fractal box"
)
raise ValueError
fractalrange = (round_value(limits[0] / grid.dy), round_value(limits[1] / grid.dy))
# yminus surface
@@ -164,7 +170,9 @@ class AddGrass(UserObjectGeometry):
elif zs == zf:
if zs not in [volume.zs, volume.zf]:
logger.exception(f"{self.__str__()} must specify external surfaces on a fractal box")
logger.exception(
f"{self.__str__()} must specify external surfaces on a fractal box"
)
raise ValueError
fractalrange = (round_value(limits[0] / grid.dz), round_value(limits[1] / grid.dz))
# zminus surface
@@ -219,7 +227,8 @@ class AddGrass(UserObjectGeometry):
# probability values, and convert the 1D index back into a x, y index
# for the original surface.
bladesindex = np.unravel_index(
np.digitize(A, probability1D), (surface.fractalsurface.shape[0], surface.fractalsurface.shape[1])
np.digitize(A, probability1D),
(surface.fractalsurface.shape[0], surface.fractalsurface.shape[1]),
)
# Set the fractal range to minimum and maximum heights of the grass blades
@@ -227,7 +236,9 @@ class AddGrass(UserObjectGeometry):
# Set the fractal surface using the pre-calculated spatial distribution
# and a random height
surface.fractalsurface = np.zeros((surface.fractalsurface.shape[0], surface.fractalsurface.shape[1]))
surface.fractalsurface = np.zeros(
(surface.fractalsurface.shape[0], surface.fractalsurface.shape[1])
)
for i in range(len(bladesindex[0])):
surface.fractalsurface[bladesindex[0][i], bladesindex[1][i]] = R.randint(
surface.fractalrange[0], surface.fractalrange[1], size=1

22
gprMax/cmds_geometry/add_surface_roughness.py
查看文件

@@ -103,7 +103,9 @@ class AddSurfaceRoughness(UserObjectGeometry):
xf, yf, zf = p2
if frac_dim < 0:
logger.exception(f"{self.__str__()} requires a positive value for the " + "fractal dimension")
logger.exception(
f"{self.__str__()} requires a positive value for the " + "fractal dimension"
)
raise ValueError
if weighting[0] < 0:
logger.exception(
@@ -124,7 +126,9 @@ class AddSurfaceRoughness(UserObjectGeometry):
logger.exception(f"{self.__str__()} dimensions are not specified correctly")
raise ValueError
if xs not in [volume.xs, volume.xf]:
logger.exception(f"{self.__str__()} can only be used on the external surfaces of a fractal box")
logger.exception(
f"{self.__str__()} can only be used on the external surfaces of a fractal box"
)
raise ValueError
fractalrange = (round_value(limits[0] / grid.dx), round_value(limits[1] / grid.dx))
# xminus surface
@@ -155,7 +159,10 @@ class AddSurfaceRoughness(UserObjectGeometry):
logger.exception(f"{self.__str__()} dimensions are not specified correctly")
raise ValueError
if ys not in [volume.ys, volume.yf]:
logger.exception(f"{self.__str__()} can only be used on the external " + "surfaces of a fractal box")
logger.exception(
f"{self.__str__()} can only be used on the external "
+ "surfaces of a fractal box"
)
raise ValueError
fractalrange = (round_value(limits[0] / grid.dy), round_value(limits[1] / grid.dy))
# yminus surface
@@ -183,7 +190,10 @@ class AddSurfaceRoughness(UserObjectGeometry):
elif zs == zf:
if zs not in [volume.zs, volume.zf]:
logger.exception(f"{self.__str__()} can only be used on the external " + "surfaces of a fractal box")
logger.exception(
f"{self.__str__()} can only be used on the external "
+ "surfaces of a fractal box"
)
raise ValueError
fractalrange = (round_value(limits[0] / grid.dz), round_value(limits[1] / grid.dz))
# zminus surface
@@ -222,7 +232,9 @@ class AddSurfaceRoughness(UserObjectGeometry):
# List of existing surfaces IDs
existingsurfaceIDs = [x.surfaceID for x in volume.fractalsurfaces]
if surface.surfaceID in existingsurfaceIDs:
logger.exception(f"{self.__str__()} has already been used on the {surface.surfaceID} surface")
logger.exception(
f"{self.__str__()} has already been used on the {surface.surfaceID} surface"
)
raise ValueError
surface.generate_fractal_surface()

17
gprMax/cmds_geometry/add_surface_water.py
查看文件

@@ -94,7 +94,9 @@ class AddSurfaceWater(UserObjectGeometry):
logger.exception(f"{self.__str__()} dimensions are not specified correctly")
raise ValueError
if xs not in [volume.xs, volume.xf]:
logger.exception(f"{self.__str__()} can only be used on the external surfaces of a fractal box")
logger.exception(
f"{self.__str__()} can only be used on the external surfaces of a fractal box"
)
raise ValueError
# xminus surface
if xs == volume.xs:
@@ -110,7 +112,9 @@ class AddSurfaceWater(UserObjectGeometry):
logger.exception(f"{self.__str__()} dimensions are not specified correctly")
raise ValueError
if ys not in [volume.ys, volume.yf]:
logger.exception(f"{self.__str__()} can only be used on the external surfaces of a fractal box")
logger.exception(
f"{self.__str__()} can only be used on the external surfaces of a fractal box"
)
raise ValueError
# yminus surface
if ys == volume.ys:
@@ -123,7 +127,9 @@ class AddSurfaceWater(UserObjectGeometry):
elif zs == zf:
if zs not in [volume.zs, volume.zf]:
logger.exception(f"{self.__str__()} can only be used on the external surfaces of a fractal box")
logger.exception(
f"{self.__str__()} can only be used on the external surfaces of a fractal box"
)
raise ValueError
# zminus surface
if zs == volume.zs:
@@ -148,7 +154,10 @@ class AddSurfaceWater(UserObjectGeometry):
surface.filldepth = filldepthcells
# Check that requested fill depth falls within range of surface roughness
if surface.filldepth < surface.fractalrange[0] or surface.filldepth > surface.fractalrange[1]:
if (
surface.filldepth < surface.fractalrange[0]
or surface.filldepth > surface.fractalrange[1]
):
logger.exception(
f"{self.__str__()} requires a value for the depth of water that lies with the "
f"range of the requested surface roughness"

2
gprMax/cmds_geometry/box.py
查看文件

@@ -23,7 +23,6 @@ import numpy as np
import gprMax.config as config
from ..cython.geometry_primitives import build_box
from ..hash_cmds_geometry import check_averaging
from ..materials import Material
from .cmds_geometry import UserObjectGeometry, rotate_2point_object
@@ -88,7 +87,6 @@ class Box(UserObjectGeometry):
try:
# Try user-specified averaging
averagebox = self.kwargs["averaging"]
averagebox = check_averaging(averagebox)
except KeyError:
# Otherwise go with the grid default
averagebox = grid.averagevolumeobjects

10
gprMax/cmds_geometry/cone.py
查看文件

@@ -21,7 +21,6 @@ import logging
import numpy as np
from ..cython.geometry_primitives import build_cone
from ..hash_cmds_geometry import check_averaging
from ..materials import Material
from .cmds_geometry import UserObjectGeometry
@@ -63,7 +62,6 @@ class Cone(UserObjectGeometry):
try:
# Try user-specified averaging
averagecone = self.kwargs["averaging"]
averagecone = check_averaging(averagecone)
except KeyError:
# Otherwise go with the grid default
averagecone = grid.averagevolumeobjects
@@ -87,11 +85,15 @@ class Cone(UserObjectGeometry):
x2, y2, z2 = uip.round_to_grid(p2)
if r1 < 0:
logger.exception(f"{self.__str__()} the radius of the first face {r1:g} should be a positive value.")
logger.exception(
f"{self.__str__()} the radius of the first face {r1:g} should be a positive value."
)
raise ValueError
if r2 < 0:
logger.exception(f"{self.__str__()} the radius of the second face {r2:g} should be a positive value.")
logger.exception(
f"{self.__str__()} the radius of the second face {r2:g} should be a positive value."
)
raise ValueError
if r1 == 0 and r2 == 0:

2
gprMax/cmds_geometry/cylinder.py
查看文件

@@ -21,7 +21,6 @@ import logging
import numpy as np
from ..cython.geometry_primitives import build_cylinder
from ..hash_cmds_geometry import check_averaging
from ..materials import Material
from .cmds_geometry import UserObjectGeometry
@@ -60,7 +59,6 @@ class Cylinder(UserObjectGeometry):
try:
# Try user-specified averaging
averagecylinder = self.kwargs["averaging"]
averagecylinder = check_averaging(averagecylinder)
except KeyError:
# Otherwise go with the grid default
averagecylinder = grid.averagevolumeobjects

10
gprMax/cmds_geometry/cylindrical_sector.py
查看文件

@@ -21,7 +21,6 @@ import logging
import numpy as np
from ..cython.geometry_primitives import build_cylindrical_sector
from ..hash_cmds_geometry import check_averaging
from ..materials import Material
from .cmds_geometry import UserObjectGeometry
@@ -75,7 +74,6 @@ class CylindricalSector(UserObjectGeometry):
try:
# Try user-specified averaging
averagecylindricalsector = self.kwargs["averaging"]
averagecylindricalsector = check_averaging(averagecylindricalsector)
except KeyError:
# Otherwise go with the grid default
averagecylindricalsector = grid.averagevolumeobjects
@@ -101,10 +99,14 @@ class CylindricalSector(UserObjectGeometry):
if r <= 0:
logger.exception(f"{self.__str__()} the radius {r:g} should be a positive value.")
if sectorstartangle < 0 or sectorangle <= 0:
logger.exception(f"{self.__str__()} the starting angle and sector angle should be a positive values.")
logger.exception(
f"{self.__str__()} the starting angle and sector angle should be a positive values."
)
raise ValueError
if sectorstartangle >= 2 * np.pi or sectorangle >= 2 * np.pi:
logger.exception(f"{self.__str__()} the starting angle and sector angle must be less than 360 degrees.")
logger.exception(
f"{self.__str__()} the starting angle and sector angle must be less than 360 degrees."
)
raise ValueError
# Look up requested materials in existing list of material instances

2
gprMax/cmds_geometry/ellipsoid.py
查看文件

@@ -21,7 +21,6 @@ import logging
import numpy as np
from ..cython.geometry_primitives import build_ellipsoid
from ..hash_cmds_geometry import check_averaging
from ..materials import Material
from .cmds_geometry import UserObjectGeometry
@@ -61,7 +60,6 @@ class Ellipsoid(UserObjectGeometry):
try:
# Try user-specified averaging
averageellipsoid = self.kwargs["averaging"]
averageellipsoid = check_averaging(averageellipsoid)
except KeyError:
# Otherwise go with the grid default
averageellipsoid = grid.averagevolumeobjects

201
gprMax/cmds_geometry/fractal_box.py
查看文件

@@ -116,16 +116,24 @@ class FractalBox(UserObjectGeometry):
xf, yf, zf = p2
if frac_dim < 0:
logger.exception(f"{self.__str__()} requires a positive value for the fractal dimension")
logger.exception(
f"{self.__str__()} requires a positive value for the fractal dimension"
)
raise ValueError
if weighting[0] < 0:
logger.exception(f"{self.__str__()} requires a positive value for the fractal weighting in the x direction")
logger.exception(
f"{self.__str__()} requires a positive value for the fractal weighting in the x direction"
)
raise ValueError
if weighting[1] < 0:
logger.exception(f"{self.__str__()} requires a positive value for the fractal weighting in the y direction")
logger.exception(
f"{self.__str__()} requires a positive value for the fractal weighting in the y direction"
)
raise ValueError
if weighting[2] < 0:
logger.exception(f"{self.__str__()} requires a positive value for the fractal weighting in the z direction")
logger.exception(
f"{self.__str__()} requires a positive value for the fractal weighting in the z direction"
)
if n_materials < 0:
logger.exception(f"{self.__str__()} requires a positive value for the number of bins")
raise ValueError
@@ -138,18 +146,25 @@ class FractalBox(UserObjectGeometry):
if mixingmodel:
if nbins == 1:
logger.exception(f"{self.__str__()} must be used with more than one material from the mixing model.")
logger.exception(
f"{self.__str__()} must be used with more than one material from the mixing model."
)
raise ValueError
if isinstance(mixingmodel, ListMaterial) and nbins > len(mixingmodel.mat):
logger.exception(
f"{self.__str__()} too many materials/bins " "requested compared to materials in " "mixing model."
f"{self.__str__()} too many materials/bins "
"requested compared to materials in "
"mixing model."
)
raise ValueError
# Create materials from mixing model as number of bins now known
# from fractal_box command.
mixingmodel.calculate_properties(nbins, grid)
elif not material:
logger.exception(f"{self.__str__()} mixing model or material with " + "ID {mixing_model_id} does not exist")
logger.exception(
f"{self.__str__()} mixing model or material with "
+ "ID {mixing_model_id} does not exist"
)
raise ValueError
self.volume = FractalVolume(xs, xf, ys, yf, zs, zf, frac_dim, seed)
@@ -221,7 +236,9 @@ class FractalBox(UserObjectGeometry):
(self.volume.nx, self.volume.ny, self.volume.nz),
dtype=config.sim_config.dtypes["float_or_double"],
)
materialnumID = next(x.numID for x in grid.materials if x.ID == self.volume.operatingonID)
materialnumID = next(
x.numID for x in grid.materials if x.ID == self.volume.operatingonID
)
self.volume.fractalvolume *= materialnumID
else:
self.volume.generate_fractal_volume()
@@ -229,7 +246,9 @@ class FractalBox(UserObjectGeometry):
for j in range(0, self.volume.ny):
for k in range(0, self.volume.nz):
numberinbin = self.volume.fractalvolume[i, j, k]
self.volume.fractalvolume[i, j, k] = self.volume.mixingmodel.matID[int(numberinbin)]
self.volume.fractalvolume[i, j, k] = self.volume.mixingmodel.matID[
int(numberinbin)
]
self.volume.generate_volume_mask()
@@ -242,28 +261,49 @@ class FractalBox(UserObjectGeometry):
for j in range(surface.ys, surface.yf):
for k in range(surface.zs, surface.zf):
if i > surface.fractalsurface[j - surface.ys, k - surface.zs]:
self.volume.mask[i - self.volume.xs, j - self.volume.ys, k - self.volume.zs] = 1
self.volume.mask[
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 1
elif surface.filldepth > 0 and i > surface.filldepth:
self.volume.mask[i - self.volume.xs, j - self.volume.ys, k - self.volume.zs] = 2
self.volume.mask[
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 2
else:
self.volume.mask[i - self.volume.xs, j - self.volume.ys, k - self.volume.zs] = 0
self.volume.mask[
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 0
elif surface.surfaceID == "xplus":
if not surface.ID:
for i in range(surface.fractalrange[0], surface.fractalrange[1]):
for j in range(surface.ys, surface.yf):
for k in range(surface.zs, surface.zf):
if i < surface.fractalsurface[j - surface.ys, k - surface.zs]:
if (
i
< surface.fractalsurface[j - surface.ys, k - surface.zs]
):
self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 1
elif surface.filldepth > 0 and i < surface.filldepth:
self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 2
else:
self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 0
elif surface.ID == "grass":
g = surface.grass[0]
@@ -275,9 +315,14 @@ class FractalBox(UserObjectGeometry):
height = 0
for i in range(self.volume.xs, surface.fractalrange[1]):
if (
i < surface.fractalsurface[j - surface.ys, k - surface.zs]
i
< surface.fractalsurface[
j - surface.ys, k - surface.zs
]
and self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
]
!= 1
):
@@ -311,11 +356,15 @@ class FractalBox(UserObjectGeometry):
i
> self.volume.originalxf
- (
surface.fractalsurface[j - surface.ys, k - surface.zs]
surface.fractalsurface[
j - surface.ys, k - surface.zs
]
- self.volume.originalxf
)
and self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
]
== 1
):
@@ -343,28 +392,49 @@ class FractalBox(UserObjectGeometry):
for j in range(surface.fractalrange[0], surface.fractalrange[1]):
for k in range(surface.zs, surface.zf):
if j > surface.fractalsurface[i - surface.xs, k - surface.zs]:
self.volume.mask[i - self.volume.xs, j - self.volume.ys, k - self.volume.zs] = 1
self.volume.mask[
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 1
elif surface.filldepth > 0 and j > surface.filldepth:
self.volume.mask[i - self.volume.xs, j - self.volume.ys, k - self.volume.zs] = 2
self.volume.mask[
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 2
else:
self.volume.mask[i - self.volume.xs, j - self.volume.ys, k - self.volume.zs] = 0
self.volume.mask[
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 0
elif surface.surfaceID == "yplus":
if not surface.ID:
for i in range(surface.xs, surface.xf):
for j in range(surface.fractalrange[0], surface.fractalrange[1]):
for k in range(surface.zs, surface.zf):
if j < surface.fractalsurface[i - surface.xs, k - surface.zs]:
if (
j
< surface.fractalsurface[i - surface.xs, k - surface.zs]
):
self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 1
elif surface.filldepth > 0 and j < surface.filldepth:
self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 2
else:
self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 0
elif surface.ID == "grass":
g = surface.grass[0]
@@ -376,9 +446,14 @@ class FractalBox(UserObjectGeometry):
height = 0
for j in range(self.volume.ys, surface.fractalrange[1]):
if (
j < surface.fractalsurface[i - surface.xs, k - surface.zs]
j
< surface.fractalsurface[
i - surface.xs, k - surface.zs
]
and self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
]
!= 1
):
@@ -412,11 +487,15 @@ class FractalBox(UserObjectGeometry):
j
> self.volume.originalyf
- (
surface.fractalsurface[i - surface.xs, k - surface.zs]
surface.fractalsurface[
i - surface.xs, k - surface.zs
]
- self.volume.originalyf
)
and self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
]
== 1
):
@@ -444,28 +523,51 @@ class FractalBox(UserObjectGeometry):
for j in range(surface.ys, surface.yf):
for k in range(surface.fractalrange[0], surface.fractalrange[1]):
if k > surface.fractalsurface[i - surface.xs, j - surface.ys]:
self.volume.mask[i - self.volume.xs, j - self.volume.ys, k - self.volume.zs] = 1
self.volume.mask[
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 1
elif surface.filldepth > 0 and k > surface.filldepth:
self.volume.mask[i - self.volume.xs, j - self.volume.ys, k - self.volume.zs] = 2
self.volume.mask[
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 2
else:
self.volume.mask[i - self.volume.xs, j - self.volume.ys, k - self.volume.zs] = 0
self.volume.mask[
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 0
elif surface.surfaceID == "zplus":
if not surface.ID:
for i in range(surface.xs, surface.xf):
for j in range(surface.ys, surface.yf):
for k in range(surface.fractalrange[0], surface.fractalrange[1]):
if k < surface.fractalsurface[i - surface.xs, j - surface.ys]:
for k in range(
surface.fractalrange[0], surface.fractalrange[1]
):
if (
k
< surface.fractalsurface[i - surface.xs, j - surface.ys]
):
self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 1
elif surface.filldepth > 0 and k < surface.filldepth:
self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 2
else:
self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
] = 0
elif surface.ID == "grass":
g = surface.grass[0]
@@ -477,9 +579,14 @@ class FractalBox(UserObjectGeometry):
height = 0
for k in range(self.volume.zs, surface.fractalrange[1]):
if (
k < surface.fractalsurface[i - surface.xs, j - surface.ys]
k
< surface.fractalsurface[
i - surface.xs, j - surface.ys
]
and self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
]
!= 1
):
@@ -513,11 +620,15 @@ class FractalBox(UserObjectGeometry):
k
> self.volume.originalzf
- (
surface.fractalsurface[i - surface.xs, j - surface.ys]
surface.fractalsurface[
i - surface.xs, j - surface.ys
]
- self.volume.originalzf
)
and self.volume.mask[
i - self.volume.xs, j - self.volume.ys, k - self.volume.zs
i - self.volume.xs,
j - self.volume.ys,
k - self.volume.zs,
]
== 1
):
@@ -575,7 +686,9 @@ class FractalBox(UserObjectGeometry):
for j in range(0, self.volume.ny):
for k in range(0, self.volume.nz):
numberinbin = self.volume.fractalvolume[i, j, k]
self.volume.fractalvolume[i, j, k] = self.volume.mixingmodel.matID[int(numberinbin)]
self.volume.fractalvolume[i, j, k] = self.volume.mixingmodel.matID[
int(numberinbin)
]
data = self.volume.fractalvolume.astype("int16", order="C")
build_voxels_from_array(

16
gprMax/cmds_geometry/geometry_objects_read.py
查看文件

@@ -126,12 +126,18 @@ class GeometryObjectsRead(UserObjectGeometry):
rigidE = f["/rigidE"][:]
rigidH = f["/rigidH"][:]
ID = f["/ID"][:]
grid.solid[xs : xs + data.shape[0], ys : ys + data.shape[1], zs : zs + data.shape[2]] = (
data + numexistmaterials
grid.solid[
xs : xs + data.shape[0], ys : ys + data.shape[1], zs : zs + data.shape[2]
] = (data + numexistmaterials)
grid.rigidE[
:, xs : xs + rigidE.shape[1], ys : ys + rigidE.shape[2], zs : zs + rigidE.shape[3]
] = rigidE
grid.rigidH[
:, xs : xs + rigidH.shape[1], ys : ys + rigidH.shape[2], zs : zs + rigidH.shape[3]
] = rigidH
grid.ID[:, xs : xs + ID.shape[1], ys : ys + ID.shape[2], zs : zs + ID.shape[3]] = (
ID + numexistmaterials
)
grid.rigidE[:, xs : xs + rigidE.shape[1], ys : ys + rigidE.shape[2], zs : zs + rigidE.shape[3]] = rigidE
grid.rigidH[:, xs : xs + rigidH.shape[1], ys : ys + rigidH.shape[2], zs : zs + rigidH.shape[3]] = rigidH
grid.ID[:, xs : xs + ID.shape[1], ys : ys + ID.shape[2], zs : zs + ID.shape[3]] = ID + numexistmaterials
logger.info(
f"{self.grid_name(grid)}Geometry objects from file {geofile} "
f"inserted at {xs * grid.dx:g}m, {ys * grid.dy:g}m, "

2
gprMax/cmds_geometry/sphere.py
查看文件

@@ -21,7 +21,6 @@ import logging
import numpy as np
from ..cython.geometry_primitives import build_sphere
from ..hash_cmds_geometry import check_averaging
from ..materials import Material
from .cmds_geometry import UserObjectGeometry
@@ -56,7 +55,6 @@ class Sphere(UserObjectGeometry):
try:
# Try user-specified averaging
averagesphere = self.kwargs["averaging"]
averagesphere = check_averaging(averagesphere)
except KeyError:
# Otherwise go with the grid default
averagesphere = grid.averagevolumeobjects

2
gprMax/cmds_geometry/triangle.py
查看文件

@@ -21,7 +21,6 @@ import logging
import numpy as np
from ..cython.geometry_primitives import build_triangle
from ..hash_cmds_geometry import check_averaging
from ..materials import Material
from .cmds_geometry import UserObjectGeometry, rotate_point
@@ -81,7 +80,6 @@ class Triangle(UserObjectGeometry):
try:
# Try user-specified averaging
averagetriangularprism = self.kwargs["averaging"]
averagetriangularprism = check_averaging(averagetriangularprism)
except KeyError:
# Otherwise go with the grid default
averagetriangularprism = grid.averagevolumeobjects

205
gprMax/cmds_multiuse.py
查看文件

@@ -26,6 +26,7 @@ from scipy import interpolate
import gprMax.config as config
from gprMax.grid.fdtd_grid import FDTDGrid
from gprMax.model import Model
from gprMax.user_inputs import MainGridUserInput
from .cmds_geometry.cmds_geometry import (
@@ -74,11 +75,11 @@ class UserObjectMulti(ABC):
return f"{self.hash}: {s[:-1]}"
@abstractmethod
def build(self, grid: FDTDGrid, uip: MainGridUserInput):
"""Creates object and adds it to grid."""
def build(self, model: Model, uip: MainGridUserInput):
"""Creates object and adds it to model."""
pass
# TODO: Check if this is actually needed
# TODO: Make _do_rotate not use a grid object
def rotate(self, axis, angle, origin=None):
"""Rotates object (specialised for each object)."""
pass
@@ -87,16 +88,19 @@ class UserObjectMulti(ABC):
"""Readable string of parameters given to object."""
return f"{self.hash}: {str(self.kwargs)}"
def grid_name(self, grid):
def grid_name(self, grid: FDTDGrid) -> str:
"""Returns subgrid name for use with logging info. Returns an empty
string if the grid is the main grid.
"""
if isinstance(grid, SubGridBaseGrid):
return f"[{grid.name}] "
else:
return ""
def model_name(self, model: Model) -> str:
"""Returns model name for use with logging info."""
return f"[{model.title}] "
class ExcitationFile(UserObjectMulti):
"""An ASCII file that contains columns of amplitude values that specify
@@ -114,7 +118,7 @@ class ExcitationFile(UserObjectMulti):
self.order = 1
self.hash = "#excitation_file"
def build(self, grid, uip):
def build(self, model, uip):
try:
kwargs = {}
excitationfile = self.kwargs["filepath"]
@@ -136,6 +140,7 @@ class ExcitationFile(UserObjectMulti):
if not excitationfile.exists():
excitationfile = Path(config.sim_config.input_file_path.parent, excitationfile)
grid = uip.grid
logger.info(self.grid_name(grid) + f"Excitation file: {excitationfile}")
# Get waveform names
@@ -153,7 +158,7 @@ class ExcitationFile(UserObjectMulti):
waveformvalues = waveformvalues[:, 1:]
timestr = "user-defined time array"
else:
waveformtime = np.arange(0, grid.timewindow + grid.dt, grid.dt)
waveformtime = np.arange(0, model.timewindow + grid.dt, grid.dt)
timestr = "simulation time array"
for i, waveformID in enumerate(waveformIDs):
@@ -213,7 +218,7 @@ class Waveform(UserObjectMulti):
self.order = 2
self.hash = "#waveform"
def build(self, grid, uip):
def build(self, model, uip):
try:
wavetype = self.kwargs["wave_type"].lower()
except KeyError:
@@ -227,6 +232,7 @@ class Waveform(UserObjectMulti):
)
raise ValueError
grid = uip.grid
if wavetype != "user":
try:
amp = self.kwargs["amp"]
@@ -234,7 +240,7 @@ class Waveform(UserObjectMulti):
ID = self.kwargs["id"]
except KeyError:
logger.exception(
self.params_str() + (" builtin waveforms " "require exactly four parameters.")
self.params_str() + (" builtin waveforms require exactly four parameters.")
)
raise
if freq <= 0:
@@ -244,7 +250,7 @@ class Waveform(UserObjectMulti):
)
raise ValueError
if any(x.ID == ID for x in grid.waveforms):
logger.exception(self.params_str() + (f" with ID {ID} already " "exists."))
logger.exception(self.params_str() + (f" with ID {ID} already exists."))
raise ValueError
w = WaveformUser()
@@ -271,14 +277,14 @@ class Waveform(UserObjectMulti):
except KeyError:
logger.exception(
self.params_str()
+ (" a user-defined " "waveform requires at least two parameters.")
+ (" a user-defined waveform requires at least two parameters.")
)
raise
if "user_time" in self.kwargs:
waveformtime = self.kwargs["user_time"]
else:
waveformtime = np.arange(0, grid.timewindow + grid.dt, grid.dt)
waveformtime = np.arange(0, model.timewindow + grid.dt, grid.dt)
# Set args for interpolation if given by user
if "kind" in self.kwargs:
@@ -287,7 +293,7 @@ class Waveform(UserObjectMulti):
kwargs["fill_value"] = self.kwargs["fill_value"]
if any(x.ID == ID for x in grid.waveforms):
logger.exception(self.params_str() + (f" with ID {ID} already " "exists."))
logger.exception(self.params_str() + (f" with ID {ID} already exists."))
raise ValueError
w = WaveformUser()
@@ -295,9 +301,7 @@ class Waveform(UserObjectMulti):
w.type = wavetype
w.userfunc = interpolate.interp1d(waveformtime, uservalues, **kwargs)
logger.info(
self.grid_name(grid) + (f"Waveform {w.ID} that is " "user-defined created.")
)
logger.info(self.grid_name(grid) + (f"Waveform {w.ID} that is user-defined created."))
grid.waveforms.append(w)
@@ -335,40 +339,41 @@ class VoltageSource(UserObjectMulti):
rot_pts = rotate_2point_object(rot_pol_pts, self.axis, self.angle, self.origin)
self.kwargs["p1"] = tuple(rot_pts[0, :])
def build(self, grid, uip):
def build(self, model, uip):
try:
p1 = self.kwargs["p1"]
polarisation = self.kwargs["polarisation"].lower()
resistance = self.kwargs["resistance"]
waveform_id = self.kwargs["waveform_id"]
except KeyError:
logger.exception(self.params_str() + (" requires at least six " "parameters."))
logger.exception(self.params_str() + (" requires at least six parameters."))
raise
grid = uip.grid
if self.do_rotate:
self._do_rotate(grid)
# Check polarity & position parameters
if polarisation not in ("x", "y", "z"):
logger.exception(self.params_str() + (" polarisation must be " "x, y, or z."))
logger.exception(self.params_str() + (" polarisation must be x, y, or z."))
raise ValueError
if "2D TMx" in config.get_model_config().mode and polarisation in [
"y",
"z",
]:
logger.exception(self.params_str() + (" polarisation must be x in " "2D TMx mode."))
logger.exception(self.params_str() + (" polarisation must be x in 2D TMx mode."))
raise ValueError
elif "2D TMy" in config.get_model_config().mode and polarisation in [
"x",
"z",
]:
logger.exception(self.params_str() + (" polarisation must be y in " "2D TMy mode."))
logger.exception(self.params_str() + (" polarisation must be y in 2D TMy mode."))
raise ValueError
elif "2D TMz" in config.get_model_config().mode and polarisation in [
"x",
"y",
]:
logger.exception(self.params_str() + (" polarisation must be z in " "2D TMz mode."))
logger.exception(self.params_str() + (" polarisation must be z in 2D TMz mode."))
raise ValueError
xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.params_str())
@@ -376,14 +381,14 @@ class VoltageSource(UserObjectMulti):
if resistance < 0:
logger.exception(
self.params_str() + (" requires a source " "resistance of zero " "or greater.")
self.params_str() + (" requires a source resistance of zero or greater.")
)
raise ValueError
# Check if there is a waveformID in the waveforms list
if not any(x.ID == waveform_id for x in grid.waveforms):
logger.exception(
self.params_str() + (" there is no waveform with " "the identifier {waveform_id}.")
self.params_str() + (" there is no waveform with the identifier {waveform_id}.")
)
raise ValueError
@@ -403,7 +408,7 @@ class VoltageSource(UserObjectMulti):
+ ")"
)
v.resistance = resistance
v.waveformID = waveform_id
v.waveform = grid.get_waveform_by_id(waveform_id)
try:
start = self.kwargs["start"]
@@ -412,37 +417,36 @@ class VoltageSource(UserObjectMulti):
if start < 0:
logger.exception(
self.params_str()
+ (" delay of the initiation " "of the source should not " "be less than zero.")
+ (" delay of the initiation of the source should not be less than zero.")
)
raise ValueError
if stop < 0:
logger.exception(
self.params_str()
+ (" time to remove the " "source should not be " "less than zero.")
self.params_str() + (" time to remove the source should not be less than zero.")
)
raise ValueError
if stop - start <= 0:
logger.exception(
self.params_str()
+ (" duration of the source " "should not be zero or " "less.")
self.params_str() + (" duration of the source should not be zero or less.")
)
raise ValueError
v.start = start
v.stop = min(stop, grid.timewindow)
v.stop = min(stop, model.timewindow)
startstop = f" start time {v.start:g} secs, finish time {v.stop:g} secs "
except KeyError:
v.start = 0
v.stop = grid.timewindow
v.stop = model.timewindow
startstop = " "
v.calculate_waveform_values(grid)
iterations = grid.iterations if isinstance(grid, SubGridBaseGrid) else model.iterations
v.calculate_waveform_values(iterations, grid.dt)
logger.info(
f"{self.grid_name(grid)}Voltage source with polarity "
f"{v.polarisation} at {p2[0]:g}m, {p2[1]:g}m, {p2[2]:g}m, "
f"resistance {v.resistance:.1f} Ohms,"
+ startstop
+ f"using waveform {v.waveformID} created."
+ f"using waveform {v.waveform.ID} created."
)
grid.voltagesources.append(v)
@@ -481,7 +485,7 @@ class HertzianDipole(UserObjectMulti):
rot_pts = rotate_2point_object(rot_pol_pts, self.axis, self.angle, self.origin)
self.kwargs["p1"] = tuple(rot_pts[0, :])
def build(self, grid, uip):
def build(self, model, uip):
try:
polarisation = self.kwargs["polarisation"].lower()
p1 = self.kwargs["p1"]
@@ -490,30 +494,31 @@ class HertzianDipole(UserObjectMulti):
logger.exception(f"{self.params_str()} requires at least 3 parameters.")
raise
grid = uip.grid
if self.do_rotate:
self._do_rotate(grid)
# Check polarity & position parameters
if polarisation not in ("x", "y", "z"):
logger.exception(self.params_str() + (" polarisation must be " "x, y, or z."))
logger.exception(self.params_str() + " polarisation must be x, y, or z.")
raise ValueError
if "2D TMx" in config.get_model_config().mode and polarisation in [
"y",
"z",
]:
logger.exception(self.params_str() + (" polarisation must be x in " "2D TMx mode."))
logger.exception(self.params_str() + " polarisation must be x in 2D TMx mode.")
raise ValueError
elif "2D TMy" in config.get_model_config().mode and polarisation in [
"x",
"z",
]:
logger.exception(self.params_str() + (" polarisation must be y in " "2D TMy mode."))
logger.exception(self.params_str() + " polarisation must be y in 2D TMy mode.")
raise ValueError
elif "2D TMz" in config.get_model_config().mode and polarisation in [
"x",
"y",
]:
logger.exception(self.params_str() + (" polarisation must be z in " "2D TMz mode."))
logger.exception(self.params_str() + " polarisation must be z in 2D TMz mode.")
raise ValueError
xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.params_str())
@@ -544,7 +549,7 @@ class HertzianDipole(UserObjectMulti):
h.ycoordorigin = ycoord
h.zcoordorigin = zcoord
h.ID = f"{h.__class__.__name__}({str(h.xcoord)},{str(h.ycoord)},{str(h.zcoord)})"
h.waveformID = waveform_id
h.waveform = grid.get_waveform_by_id(waveform_id)
try:
# Check source start & source remove time parameters
@@ -566,14 +571,15 @@ class HertzianDipole(UserObjectMulti):
)
raise ValueError
h.start = start
h.stop = min(stop, grid.timewindow)
h.stop = min(stop, model.timewindow)
startstop = f" start time {h.start:g} secs, finish time {h.stop:g} secs "
except KeyError:
h.start = 0
h.stop = grid.timewindow
h.stop = model.timewindow
startstop = " "
h.calculate_waveform_values(grid)
iterations = grid.iterations if isinstance(grid, SubGridBaseGrid) else model.iterations
h.calculate_waveform_values(iterations, grid.dt)
if config.get_model_config().mode == "2D":
logger.info(
@@ -581,7 +587,7 @@ class HertzianDipole(UserObjectMulti):
+ f"source in 2D with polarity {h.polarisation} at "
+ f"{p2[0]:g}m, {p2[1]:g}m, {p2[2]:g}m,"
+ startstop
+ f"using waveform {h.waveformID} created."
+ f"using waveform {h.waveform.ID} created."
)
else:
logger.info(
@@ -589,7 +595,7 @@ class HertzianDipole(UserObjectMulti):
+ f"polarity {h.polarisation} at {p2[0]:g}m, "
+ f"{p2[1]:g}m, {p2[2]:g}m,"
+ startstop
+ f"using waveform {h.waveformID} created."
+ f"using waveform {h.waveform.ID} created."
)
grid.hertziandipoles.append(h)
@@ -628,7 +634,7 @@ class MagneticDipole(UserObjectMulti):
rot_pts = rotate_2point_object(rot_pol_pts, self.axis, self.angle, self.origin)
self.kwargs["p1"] = tuple(rot_pts[0, :])
def build(self, grid, uip):
def build(self, model, uip):
try:
polarisation = self.kwargs["polarisation"].lower()
p1 = self.kwargs["p1"]
@@ -637,30 +643,31 @@ class MagneticDipole(UserObjectMulti):
logger.exception(f"{self.params_str()} requires at least five parameters.")
raise
grid = uip.grid
if self.do_rotate:
self._do_rotate(grid)
# Check polarity & position parameters
if polarisation not in ("x", "y", "z"):
logger.exception(self.params_str() + (" polarisation must be " "x, y, or z."))
logger.exception(self.params_str() + " polarisation must be x, y, or z.")
raise ValueError
if "2D TMx" in config.get_model_config().mode and polarisation in [
"y",
"z",
]:
logger.exception(self.params_str() + (" polarisation must be x in " "2D TMx mode."))
logger.exception(self.params_str() + " polarisation must be x in 2D TMx mode.")
raise ValueError
elif "2D TMy" in config.get_model_config().mode and polarisation in [
"x",
"z",
]:
logger.exception(self.params_str() + (" polarisation must be y in " "2D TMy mode."))
logger.exception(self.params_str() + " polarisation must be y in 2D TMy mode.")
raise ValueError
elif "2D TMz" in config.get_model_config().mode and polarisation in [
"x",
"y",
]:
logger.exception(self.params_str() + (" polarisation must be z in " "2D TMz mode."))
logger.exception(self.params_str() + " polarisation must be z in 2D TMz mode.")
raise ValueError
xcoord, ycoord, zcoord = uip.check_src_rx_point(p1, self.params_str())
@@ -691,7 +698,7 @@ class MagneticDipole(UserObjectMulti):
+ str(m.zcoord)
+ ")"
)
m.waveformID = waveform_id
m.waveform = grid.get_waveform_by_id(waveform_id)
try:
# Check source start & source remove time parameters
@@ -700,36 +707,35 @@ class MagneticDipole(UserObjectMulti):
if start < 0:
logger.exception(
self.params_str()
+ (" delay of the initiation " "of the source should not " "be less than zero.")
+ " delay of the initiation of the source should not be less than zero."
)
raise ValueError
if stop < 0:
logger.exception(
self.params_str()
+ (" time to remove the " "source should not be " "less than zero.")
self.params_str() + " time to remove the source should not be less than zero."
)
raise ValueError
if stop - start <= 0:
logger.exception(
self.params_str()
+ (" duration of the source " "should not be zero or " "less.")
self.params_str() + " duration of the source should not be zero or less."
)
raise ValueError
m.start = start
m.stop = min(stop, grid.timewindow)
m.stop = min(stop, model.timewindow)
startstop = f" start time {m.start:g} secs, finish time {m.stop:g} secs "
except KeyError:
m.start = 0
m.stop = grid.timewindow
m.stop = model.timewindow
startstop = " "
m.calculate_waveform_values(grid)
iterations = grid.iterations if isinstance(grid, SubGridBaseGrid) else model.iterations
m.calculate_waveform_values(iterations, grid.dt)
logger.info(
f"{self.grid_name(grid)}Magnetic dipole with polarity "
f"{m.polarisation} at {p2[0]:g}m, {p2[1]:g}m, {p2[2]:g}m,"
+ startstop
+ f"using waveform {m.waveformID} created."
+ f"using waveform {m.waveform.ID} created."
)
grid.magneticdipoles.append(m)
@@ -769,7 +775,7 @@ class TransmissionLine(UserObjectMulti):
rot_pts = rotate_2point_object(rot_pol_pts, self.axis, self.angle, self.origin)
self.kwargs["p1"] = tuple(rot_pts[0, :])
def build(self, grid, uip):
def build(self, model, uip):
try:
polarisation = self.kwargs["polarisation"].lower()
p1 = self.kwargs["p1"]
@@ -779,6 +785,7 @@ class TransmissionLine(UserObjectMulti):
logger.exception(f"{self.params_str()} requires at least six parameters.")
raise
grid = uip.grid
if self.do_rotate:
self._do_rotate(grid)
@@ -832,7 +839,8 @@ class TransmissionLine(UserObjectMulti):
)
raise ValueError
t = TransmissionLineUser(grid)
iterations = grid.iterations if isinstance(grid, SubGridBaseGrid) else model.iterations
t = TransmissionLineUser(iterations, grid.dt)
t.polarisation = polarisation
t.xcoord = xcoord
t.ycoord = ycoord
@@ -848,7 +856,7 @@ class TransmissionLine(UserObjectMulti):
+ ")"
)
t.resistance = resistance
t.waveformID = waveform_id
t.waveform = grid.get_waveform_by_id(waveform_id)
try:
# Check source start & source remove time parameters
@@ -873,14 +881,14 @@ class TransmissionLine(UserObjectMulti):
)
raise ValueError
t.start = start
t.stop = min(stop, grid.timewindow)
t.stop = min(stop, model.timewindow)
startstop = f" start time {t.start:g} secs, finish time {t.stop:g} secs "
except KeyError:
t.start = 0
t.stop = grid.timewindow
t.stop = model.timewindow
startstop = " "
t.calculate_waveform_values(grid)
t.calculate_waveform_values(iterations, grid.dt)
t.calculate_incident_V_I(grid)
logger.info(
@@ -888,7 +896,7 @@ class TransmissionLine(UserObjectMulti):
+ f"{t.polarisation} at {p2[0]:g}m, {p2[1]:g}m, "
+ f"{p2[2]:g}m, resistance {t.resistance:.1f} Ohms,"
+ startstop
+ f"using waveform {t.waveformID} created."
+ f"using waveform {t.waveform.ID} created."
)
grid.transmissionlines.append(t)
@@ -941,13 +949,14 @@ class Rx(UserObjectMulti):
except KeyError:
pass
def build(self, grid, uip):
def build(self, model, uip):
try:
p1 = self.kwargs["p1"]
except KeyError:
logger.exception(self.params_str())
raise
grid = uip.grid
if self.do_rotate:
self._do_rotate(grid)
@@ -958,6 +967,8 @@ class Rx(UserObjectMulti):
r.xcoord, r.ycoord, r.zcoord = p
r.xcoordorigin, r.ycoordorigin, r.zcoordorigin = p
iterations = grid.iterations if isinstance(grid, SubGridBaseGrid) else model.iterations
try:
r.ID = self.kwargs["id"]
outputs = self.kwargs["outputs"]
@@ -966,7 +977,7 @@ class Rx(UserObjectMulti):
r.ID = f"{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=config.sim_config.dtypes["float_or_double"]
iterations, dtype=config.sim_config.dtypes["float_or_double"]
)
else:
outputs.sort()
@@ -979,7 +990,7 @@ class Rx(UserObjectMulti):
for field in outputs:
if field in allowableoutputs:
r.outputs[field] = np.zeros(
grid.iterations, dtype=config.sim_config.dtypes["float_or_double"]
iterations, dtype=config.sim_config.dtypes["float_or_double"]
)
else:
logger.exception(
@@ -1015,7 +1026,7 @@ class RxArray(UserObjectMulti):
self.order = 8
self.hash = "#rx_array"
def build(self, grid, uip):
def build(self, model, uip):
try:
p1 = self.kwargs["p1"]
p2 = self.kwargs["p2"]
@@ -1063,6 +1074,7 @@ class RxArray(UserObjectMulti):
)
raise ValueError
grid = uip.grid
logger.info(
f"{self.grid_name(grid)}Receiver array "
f"{p3[0]:g}m, {p3[1]:g}m, {p3[2]:g}m, to "
@@ -1070,6 +1082,8 @@ class RxArray(UserObjectMulti):
f"{dx * grid.dx:g}m, {dy * grid.dy:g}m, {dz * grid.dz:g}m"
)
iterations = grid.iterations if isinstance(grid, SubGridBaseGrid) else model.iterations
for x in range(xs, xf + 1, dx):
for y in range(ys, yf + 1, dy):
for z in range(zs, zf + 1, dz):
@@ -1087,7 +1101,7 @@ class RxArray(UserObjectMulti):
r.ID = f"{r.__class__.__name__}({str(x)},{str(y)},{str(z)})"
for key in RxUser.defaultoutputs:
r.outputs[key] = np.zeros(
grid.iterations, dtype=config.sim_config.dtypes["float_or_double"]
iterations, dtype=config.sim_config.dtypes["float_or_double"]
)
logger.info(
f" Receiver at {p5[0]:g}m, {p5[1]:g}m, "
@@ -1123,7 +1137,8 @@ class Snapshot(UserObjectMulti):
self.order = 9
self.hash = "#snapshot"
def build(self, grid, uip):
def build(self, model, uip):
grid = uip.grid
if isinstance(grid, SubGridBaseGrid):
logger.exception(f"{self.params_str()} do not add snapshots to subgrids.")
raise ValueError
@@ -1203,7 +1218,7 @@ class Snapshot(UserObjectMulti):
f"{self.params_str()} the step size should not be less than the spatial discretisation."
)
raise ValueError
if iterations <= 0 or iterations > grid.iterations:
if iterations <= 0 or iterations > model.iterations:
logger.exception(f"{self.params_str()} time value is not valid.")
raise ValueError
@@ -1252,7 +1267,7 @@ class Material(UserObjectMulti):
self.order = 10
self.hash = "#material"
def build(self, grid, uip):
def build(self, model, uip):
try:
er = self.kwargs["er"]
se = self.kwargs["se"]
@@ -1285,6 +1300,8 @@ class Material(UserObjectMulti):
if sm < 0:
logger.exception(f"{self.params_str()} requires a positive value for magnetic loss.")
raise ValueError
grid = uip.grid
if any(x.ID == material_id for x in grid.materials):
logger.exception(f"{self.params_str()} with ID {material_id} already exists")
raise ValueError
@@ -1302,7 +1319,7 @@ class Material(UserObjectMulti):
m.er = er
logger.info(
f"{self.grid_name(grid)}Material {m.ID} with eps_r={m.er:g}, "
f"{self.model_name(model)}Material {m.ID} with eps_r={m.er:g}, "
f"sigma={m.se:g} S/m; mu_r={m.mr:g}, sigma*={m.sm:g} Ohm/m "
f"created."
)
@@ -1329,7 +1346,7 @@ class AddDebyeDispersion(UserObjectMulti):
self.order = 11
self.hash = "#add_dispersion_debye"
def build(self, grid, uip):
def build(self, model, uip):
try:
poles = self.kwargs["poles"]
er_delta = self.kwargs["er_delta"]
@@ -1344,6 +1361,7 @@ class AddDebyeDispersion(UserObjectMulti):
raise ValueError
# Look up requested materials in existing list of material instances
grid = uip.grid
materials = [y for x in material_ids for y in grid.materials if y.ID == x]
if len(materials) != len(material_ids):
@@ -1405,7 +1423,7 @@ class AddLorentzDispersion(UserObjectMulti):
self.order = 12
self.hash = "#add_dispersion_lorentz"
def build(self, grid, uip):
def build(self, model, uip):
try:
poles = self.kwargs["poles"]
er_delta = self.kwargs["er_delta"]
@@ -1421,6 +1439,7 @@ class AddLorentzDispersion(UserObjectMulti):
raise ValueError
# Look up requested materials in existing list of material instances
grid = uip.grid
materials = [y for x in material_ids for y in grid.materials if y.ID == x]
if len(materials) != len(material_ids):
@@ -1484,7 +1503,7 @@ class AddDrudeDispersion(UserObjectMulti):
self.order = 13
self.hash = "#add_dispersion_drude"
def build(self, grid, uip):
def build(self, model, uip):
try:
poles = self.kwargs["poles"]
omega = self.kwargs["omega"]
@@ -1499,6 +1518,7 @@ class AddDrudeDispersion(UserObjectMulti):
raise ValueError
# Look up requested materials in existing list of material instances
grid = uip.grid
materials = [y for x in material_ids for y in grid.materials if y.ID == x]
if len(materials) != len(material_ids):
@@ -1563,7 +1583,7 @@ class SoilPeplinski(UserObjectMulti):
self.order = 14
self.hash = "#soil_peplinski"
def build(self, grid, uip):
def build(self, model, uip):
try:
sand_fraction = self.kwargs["sand_fraction"]
clay_fraction = self.kwargs["clay_fraction"]
@@ -1606,6 +1626,7 @@ class SoilPeplinski(UserObjectMulti):
"fraction."
)
raise ValueError
grid = uip.grid
if any(x.ID == ID for x in grid.mixingmodels):
logger.exception(f"{self.params_str()} with ID {ID} already exists")
raise ValueError
@@ -1652,7 +1673,7 @@ class MaterialRange(UserObjectMulti):
self.order = 15
self.hash = "#material_range"
def build(self, grid, uip):
def build(self, model, uip):
try:
er_lower = self.kwargs["er_lower"]
er_upper = self.kwargs["er_upper"]
@@ -1710,7 +1731,7 @@ class MaterialRange(UserObjectMulti):
logger.exception(
f"{self.params_str()} requires a positive value for the upper range of magnetic loss."
)
grid = uip.grid
if any(x.ID == ID for x in grid.mixingmodels):
logger.exception(f"{self.params_str()} with ID {ID} already exists")
raise ValueError
@@ -1746,14 +1767,14 @@ class MaterialList(UserObjectMulti):
self.order = 16
self.hash = "#material_list"
def build(self, grid, uip):
def build(self, model, uip):
try:
list_of_materials = self.kwargs["list_of_materials"]
ID = self.kwargs["id"]
except KeyError:
logger.exception(f"{self.params_str()} requires at at least 2 parameters.")
raise
grid = uip.grid
if any(x.ID == ID for x in grid.mixingmodels):
logger.exception(f"{self.params_str()} with ID {ID} already exists")
raise ValueError
@@ -1801,7 +1822,7 @@ class GeometryView(UserObjectMulti):
return GeometryViewUser
def build(self, grid, uip):
def build(self, model, uip):
try:
p1 = self.kwargs["p1"]
p2 = self.kwargs["p2"]
@@ -1824,6 +1845,7 @@ class GeometryView(UserObjectMulti):
xs, ys, zs = p1
xf, yf, zf = p2
grid = uip.grid
dx, dy, dz = uip.discretise_static_point(dl)
if dx < 0 or dy < 0 or dz < 0:
@@ -1865,7 +1887,7 @@ class GeometryView(UserObjectMulti):
f"base {g.filename} created."
)
grid.geometryviews.append(g)
model.geometryviews.append(g)
class GeometryObjectsWrite(UserObjectMulti):
@@ -1886,7 +1908,11 @@ class GeometryObjectsWrite(UserObjectMulti):
self.order = 18
self.hash = "#geometry_objects_write"
def build(self, grid, uip):
def build(self, model, uip):
grid = uip.grid
if isinstance(grid, SubGridBaseGrid):
logger.exception(f"{self.params_str()} do not add geometry objects to subgrids.")
raise ValueError
try:
p1 = self.kwargs["p1"]
p2 = self.kwargs["p2"]
@@ -1910,7 +1936,7 @@ class GeometryObjectsWrite(UserObjectMulti):
f"{g.filename_materials}"
)
grid.geometryobjectswrite.append(g)
model.geometryobjects.append(g)
class PMLCFS(UserObjectMulti):
@@ -1942,7 +1968,7 @@ class PMLCFS(UserObjectMulti):
super().__init__(**kwargs)
self.order = 19
def build(self, grid, uip):
def build(self, model, uip):
try:
alphascalingprofile = self.kwargs["alphascalingprofile"]
alphascalingdirection = self.kwargs["alphascalingdirection"]
@@ -2028,6 +2054,7 @@ class PMLCFS(UserObjectMulti):
f"{cfssigma.min:g}, max: {cfssigma.max}) created."
)
grid = uip.grid
grid.pmls["cfs"].append(cfs)
if len(grid.pmls["cfs"]) > 2:

88
gprMax/cmds_singleuse.py
查看文件

@@ -15,7 +15,6 @@
#
# You should have received a copy of the GNU General Public License
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
import logging
from abc import ABC, abstractmethod
@@ -23,6 +22,7 @@ import numpy as np
import gprMax.config as config
from gprMax.grid.fdtd_grid import FDTDGrid
from gprMax.model import Model
from gprMax.user_inputs import MainGridUserInput
from .pml import PML
@@ -51,7 +51,7 @@ class UserObjectSingle(ABC):
setattr(self.props, k, v)
@abstractmethod
def build(self, grid: FDTDGrid, uip: MainGridUserInput):
def build(self, model: Model, uip: MainGridUserInput):
pass
# TODO: Check if this is actually needed
@@ -70,11 +70,11 @@ class Title(UserObjectSingle):
super().__init__(**kwargs)
self.order = 1
def build(self, G, uip):
def build(self, model, uip):
try:
title = self.kwargs["name"]
G.title = title
logger.info(f"Model title: {G.title}")
model.title = title
logger.info(f"Model title: {model.title}")
except KeyError:
pass
@@ -90,34 +90,33 @@ class Discretisation(UserObjectSingle):
super().__init__(**kwargs)
self.order = 2
def build(self, G, uip):
def build(self, model, uip):
try:
G.dl = np.array(self.kwargs["p1"])
G.dx, G.dy, G.dz = self.kwargs["p1"]
model.dl = np.array(self.kwargs["p1"], dtype=float)
except KeyError:
logger.exception(f"{self.__str__()} discretisation requires a point")
raise
if G.dl[0] <= 0:
if model.dl[0] <= 0:
logger.exception(
f"{self.__str__()} discretisation requires the "
f"x-direction spatial step to be greater than zero"
)
raise ValueError
if G.dl[1] <= 0:
if model.dl[1] <= 0:
logger.exception(
f"{self.__str__()} discretisation requires the "
f"y-direction spatial step to be greater than zero"
)
raise ValueError
if G.dl[2] <= 0:
if model.dl[2] <= 0:
logger.exception(
f"{self.__str__()} discretisation requires the "
f"z-direction spatial step to be greater than zero"
)
raise ValueError
logger.info(f"Spatial discretisation: {G.dl[0]:g} x {G.dl[1]:g} x {G.dl[2]:g}m")
logger.info(f"Spatial discretisation: {model.dl[0]:g} x {model.dl[1]:g} x {model.dl[2]:g}m")
class Domain(UserObjectSingle):
@@ -131,33 +130,37 @@ class Domain(UserObjectSingle):
super().__init__(**kwargs)
self.order = 3
def build(self, G, uip):
def build(self, model, uip):
try:
G.nx, G.ny, G.nz = uip.discretise_point(self.kwargs["p1"])
model.gnx, model.gny, model.gnz = uip.discretise_point(self.kwargs["p1"])
model.G.nx = model.gnx
model.G.ny = model.gny
model.G.nz = model.gnz
except KeyError:
logger.exception(f"{self.__str__()} please specify a point")
raise
if G.nx == 0 or G.ny == 0 or G.nz == 0:
logger.exception(f"{self.__str__()} requires at least one cell in " f"every dimension")
if model.gnx == 0 or model.gny == 0 or model.gnz == 0:
logger.exception(f"{self.__str__()} requires at least one cell in every dimension")
raise ValueError
logger.info(
f"Domain size: {self.kwargs['p1'][0]:g} x {self.kwargs['p1'][1]:g} x "
+ f"{self.kwargs['p1'][2]:g}m ({G.nx:d} x {G.ny:d} x {G.nz:d} = "
+ f"{(G.nx * G.ny * G.nz):g} cells)"
+ f"{self.kwargs['p1'][2]:g}m ({model.gnx:d} x {model.gny:d} x {model.gnz:d} = "
+ f"{(model.gnx * model.gny * model.gnz):g} cells)"
)
# Calculate time step at CFL limit; switch off appropriate PMLs for 2D
if G.nx == 1:
G = model.G
if model.gnx == 1:
config.get_model_config().mode = "2D TMx"
G.pmls["thickness"]["x0"] = 0
G.pmls["thickness"]["xmax"] = 0
elif G.ny == 1:
elif model.gny == 1:
config.get_model_config().mode = "2D TMy"
G.pmls["thickness"]["y0"] = 0
G.pmls["thickness"]["ymax"] = 0
elif G.nz == 1:
elif model.gnz == 1:
config.get_model_config().mode = "2D TMz"
G.pmls["thickness"]["z0"] = 0
G.pmls["thickness"]["zmax"] = 0
@@ -187,7 +190,7 @@ class TimeStepStabilityFactor(UserObjectSingle):
super().__init__(**kwargs)
self.order = 4
def build(self, G, uip):
def build(self, model, uip):
try:
f = self.kwargs["f"]
except KeyError:
@@ -201,10 +204,10 @@ class TimeStepStabilityFactor(UserObjectSingle):
)
raise ValueError
G.dt_mod = f
G.dt = G.dt * G.dt_mod
model.dt_mod = f
model.dt *= model.dt_mod
logger.info(f"Time step (modified): {G.dt:g} secs")
logger.info(f"Time step (modified): {model.dt:g} secs")
class TimeWindow(UserObjectSingle):
@@ -219,33 +222,33 @@ class TimeWindow(UserObjectSingle):
super().__init__(**kwargs)
self.order = 5
def build(self, G, uip):
def build(self, model, uip):
# If number of iterations given
# The +/- 1 used in calculating the number of iterations is to account for
# the fact that the solver (iterations) loop runs from 0 to < G.iterations
try:
iterations = int(self.kwargs["iterations"])
G.timewindow = (iterations - 1) * G.dt
G.iterations = iterations
model.timewindow = (iterations - 1) * model.dt
model.iterations = iterations
except KeyError:
pass
try:
tmp = float(self.kwargs["time"])
if tmp > 0:
G.timewindow = tmp
G.iterations = int(np.ceil(tmp / G.dt)) + 1
model.timewindow = tmp
model.iterations = int(np.ceil(tmp / model.dt)) + 1
else:
logger.exception(self.__str__() + " must have a value greater than zero")
raise ValueError
except KeyError:
pass
if not G.timewindow:
if not model.timewindow:
logger.exception(self.__str__() + " specify a time or number of iterations")
raise ValueError
logger.info(f"Time window: {G.timewindow:g} secs ({G.iterations} iterations)")
logger.info(f"Time window: {model.timewindow:g} secs ({model.iterations} iterations)")
class OMPThreads(UserObjectSingle):
@@ -260,7 +263,7 @@ class OMPThreads(UserObjectSingle):
super().__init__(**kwargs)
self.order = 6
def build(self, G, uip):
def build(self, model, uip):
try:
n = self.kwargs["n"]
except KeyError:
@@ -296,7 +299,8 @@ class PMLProps(UserObjectSingle):
super().__init__(**kwargs)
self.order = 7
def build(self, G, uip):
def build(self, model, uip):
G = model.G
try:
G.pmls["formulation"] = self.kwargs["formulation"]
if G.pmls["formulation"] not in PML.formulations:
@@ -348,16 +352,16 @@ class SrcSteps(UserObjectSingle):
super().__init__(**kwargs)
self.order = 8
def build(self, G, uip):
def build(self, model, uip):
try:
G.srcsteps = uip.discretise_point(self.kwargs["p1"])
model.srcsteps = uip.discretise_point(self.kwargs["p1"])
except KeyError:
logger.exception(f"{self.__str__()} requires exactly three parameters")
raise
logger.info(
f"Simple sources will step {G.srcsteps[0] * G.dx:g}m, "
f"{G.srcsteps[1] * G.dy:g}m, {G.srcsteps[2] * G.dz:g}m "
f"Simple sources will step {model.srcsteps[0] * model.dx:g}m, "
f"{model.srcsteps[1] * model.dy:g}m, {model.srcsteps[2] * model.dz:g}m "
"for each model run."
)
@@ -373,16 +377,16 @@ class RxSteps(UserObjectSingle):
super().__init__(**kwargs)
self.order = 9
def build(self, G, uip):
def build(self, model, uip):
try:
G.rxsteps = uip.discretise_point(self.kwargs["p1"])
model.rxsteps = uip.discretise_point(self.kwargs["p1"])
except KeyError:
logger.exception(f"{self.__str__()} requires exactly three parameters")
raise
logger.info(
f"All receivers will step {G.rxsteps[0] * G.dx:g}m, "
f"{G.rxsteps[1] * G.dy:g}m, {G.rxsteps[2] * G.dz:g}m "
f"All receivers will step {model.rxsteps[0] * model.dx:g}m, "
f"{model.rxsteps[1] * model.dy:g}m, {model.rxsteps[2] * model.dz:g}m "
"for each model run."
)

4
gprMax/contexts.py
查看文件

@@ -106,12 +106,12 @@ class Context:
if not model_config.reuse_geometry():
scene = self._get_scene(model_num)
model = self._create_model()
scene.create_internal_objects(model.G)
scene.create_internal_objects(model)
model.build()
if not config.sim_config.geometry_only:
solver = create_solver(model.G)
solver = create_solver(model)
model.solve(solver)
del solver

85
gprMax/fields_outputs.py
查看文件

@@ -17,15 +17,18 @@
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
import logging
from pathlib import Path
import h5py
from gprMax.grid.fdtd_grid import FDTDGrid
from ._version import __version__
logger = logging.getLogger(__name__)
def store_outputs(G):
def store_outputs(G: FDTDGrid, iteration: int):
"""Stores field component values for every receiver and transmission line.
Args:
@@ -33,7 +36,6 @@ def store_outputs(G):
"""
# Assign iteration and fields to local variables
iteration = G.iteration
Ex, Ey, Ez, Hx, Hy, Hz = G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz
for rx in G.rxs:
@@ -45,6 +47,7 @@ def store_outputs(G):
# Store current component
else:
func = globals()[output]
# TODO: Test this still works - likely broken by removing 'Ix' etc functions
rx.outputs[output][iteration] = func(rx.xcoord, rx.ycoord, rx.zcoord, Hx, Hy, Hz, G)
for tl in G.transmissionlines:
@@ -52,26 +55,29 @@ def store_outputs(G):
tl.Itotal[iteration] = tl.current[tl.antpos]
def write_hdf5_outputfile(outputfile, G):
# TODO: Add type information for grid (without a circular dependency)
def write_hdf5_outputfile(outputfile: Path, title: str, model):
"""Writes an output file in HDF5 (.h5) format.
Args:
outputfile: string of the name of the output file.
G: FDTDGrid class describing a grid in a model.
"""
# Create output file and write top-level meta data, meta data for main grid,
# and any outputs in the main grid
f = h5py.File(outputfile, "w")
f.attrs["gprMax"] = __version__
f.attrs["Title"] = G.title
write_hd5_data(f, G)
f.attrs["Title"] = title
f.attrs["Iterations"] = model.iterations
f.attrs["srcsteps"] = model.srcsteps
f.attrs["rxsteps"] = model.rxsteps
write_hd5_data(f, model.G)
# Write meta data and data for any subgrids
sg_rxs = [True for sg in G.subgrids if sg.rxs]
sg_tls = [True for sg in G.subgrids if sg.transmissionlines]
sg_rxs = [True for sg in model.subgrids if sg.rxs]
sg_tls = [True for sg in model.subgrids if sg.transmissionlines]
if sg_rxs or sg_tls:
for sg in G.subgrids:
for sg in model.subgrids:
grp = f.create_group(f"/subgrids/{sg.name}")
write_hd5_data(grp, sg, is_subgrid=True)
@@ -88,18 +94,20 @@ def write_hd5_data(basegrp, grid, is_subgrid=False):
"""
# Write meta data for grid
basegrp.attrs["Iterations"] = grid.iterations
basegrp.attrs["nx_ny_nz"] = (grid.nx, grid.ny, grid.nz)
basegrp.attrs["dx_dy_dz"] = (grid.dx, grid.dy, grid.dz)
basegrp.attrs["dt"] = grid.dt
nsrc = len(grid.voltagesources + grid.hertziandipoles + grid.magneticdipoles + grid.transmissionlines)
nsrc = len(
grid.voltagesources + grid.hertziandipoles + grid.magneticdipoles + grid.transmissionlines
)
basegrp.attrs["nsrc"] = nsrc
basegrp.attrs["nrx"] = len(grid.rxs)
basegrp.attrs["srcsteps"] = grid.srcsteps
basegrp.attrs["rxsteps"] = grid.rxsteps
if is_subgrid:
# Write additional meta data about subgrid
basegrp.attrs["Iterations"] = grid.iterations
basegrp.attrs["srcsteps"] = grid.srcsteps
basegrp.attrs["rxsteps"] = grid.rxsteps
basegrp.attrs["is_os_sep"] = grid.is_os_sep
basegrp.attrs["pml_separation"] = grid.pml_separation
basegrp.attrs["subgrid_pml_thickness"] = grid.pmls["thickness"]["x0"]
@@ -137,54 +145,3 @@ def write_hd5_data(basegrp, grid, is_subgrid=False):
for output in rx.outputs:
basegrp["rxs/rx" + str(rxindex + 1) + "/" + output] = rx.outputs[output]
def Ix(x, y, z, Hx, Hy, Hz, G):
"""Calculates the x-component of current at a grid position.
Args:
x, y, z: floats for coordinates of position in grid.
Hx, Hy, Hz: numpy array of magnetic field values.
G: FDTDGrid class describing a grid in a model.
"""
if y == 0 or z == 0:
Ix = 0
else:
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, Hy, Hz, G):
"""Calculates the y-component of current at a grid position.
Args:
x, y, z: floats for coordinates of position in grid.
Hx, Hy, Hz: numpy array of magnetic field values.
G: FDTDGrid class describing a grid in a model.
"""
if x == 0 or z == 0:
Iy = 0
else:
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, Hz, G):
"""Calculates the z-component of current at a grid position.
Args:
x, y, z: floats for coordinates of position in grid.
Hx, Hy, Hz: numpy array of magnetic field values.
G: FDTDGrid class describing a grid in a model.
"""
if x == 0 or y == 0:
Iz = 0
else:
Iz = G.dx * (Hx[x, y - 1, z] - Hx[x, y, z]) + G.dy * (Hy[x, y, z] - Hy[x - 1, y, z])
return Iz

30
gprMax/fractals.py
查看文件

@@ -56,7 +56,9 @@ class FractalSurface:
self.nz = zf - zs
self.dtype = np.dtype(np.complex128)
self.seed = seed
self.dimension = dimension # Fractal dimension from: http://dx.doi.org/10.1017/CBO9781139174695
self.dimension = (
dimension # Fractal dimension from: http://dx.doi.org/10.1017/CBO9781139174695
)
self.weighting = np.array([1, 1], dtype=np.float64)
self.fractalrange = (0, 0)
self.filldepth = 0
@@ -82,7 +84,9 @@ class FractalSurface:
self.fractalsurface = np.zeros(surfacedims, dtype=self.dtype)
# 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
rng = np.random.default_rng(seed=self.seed)
@@ -157,7 +161,9 @@ class FractalVolume:
self.averaging = False
self.dtype = np.dtype(np.complex128)
self.seed = seed
self.dimension = dimension # Fractal dimension from: http://dx.doi.org/10.1017/CBO9781139174695
self.dimension = (
dimension # Fractal dimension from: http://dx.doi.org/10.1017/CBO9781139174695
)
self.weighting = np.array([1, 1, 1], dtype=np.float64)
self.nbins = 0
self.fractalsurfaces = []
@@ -176,7 +182,9 @@ class FractalVolume:
filterscaling = np.amin(np.array([self.nx, self.ny])) / np.array([self.nx, self.ny])
filterscaling = np.insert(filterscaling, 2, 1)
else:
filterscaling = np.amin(np.array([self.nx, self.ny, self.nz])) / np.array([self.nx, self.ny, self.nz])
filterscaling = np.amin(np.array([self.nx, self.ny, self.nz])) / np.array(
[self.nx, self.ny, self.nz]
)
# Adjust weighting to account for filter scaling
self.weighting = np.multiply(self.weighting, filterscaling)
@@ -185,7 +193,11 @@ class FractalVolume:
# Positional vector at centre of array, scaled by weighting
v1 = np.array(
[self.weighting[0] * self.nx / 2, self.weighting[1] * self.ny / 2, self.weighting[2] * self.nz / 2]
[
self.weighting[0] * self.nx / 2,
self.weighting[1] * self.ny / 2,
self.weighting[2] * self.nz / 2,
]
)
# 3D array of random numbers to be convolved with the fractal function
@@ -225,7 +237,9 @@ class FractalVolume:
bins = np.linspace(np.amin(self.fractalvolume), np.amax(self.fractalvolume), self.nbins)
for j in range(self.ny):
for k in range(self.nz):
self.fractalvolume[:, j, k] = np.digitize(self.fractalvolume[:, j, k], bins, right=True)
self.fractalvolume[:, j, k] = np.digitize(
self.fractalvolume[:, j, k], bins, right=True
)
def generate_volume_mask(self):
"""Generate a 3D volume to use as a mask for adding rough surfaces,
@@ -254,7 +268,9 @@ class Grass:
"""
self.numblades = numblades
self.geometryparams = np.zeros((self.numblades, 6), dtype=config.sim_config.dtypes["float_or_double"])
self.geometryparams = np.zeros(
(self.numblades, 6), dtype=config.sim_config.dtypes["float_or_double"]
)
self.seed = seed
self.set_geometry_parameters()

81
gprMax/geometry_outputs.py
查看文件

@@ -28,6 +28,7 @@ from evtk.vtk import VtkGroup, VtkImageData, VtkUnstructuredGrid
from tqdm import tqdm
import gprMax.config as config
from gprMax.grid.fdtd_grid import FDTDGrid
from ._version import __version__
from .cython.geometry_outputs import write_lines
@@ -52,7 +53,8 @@ def save_geometry_views(gvs):
total=gv.nbytes,
unit="byte",
unit_scale=True,
desc=f"Writing geometry view file {i + 1}/{len(gvs)}, " f"{gv.filename.name}{gv.vtkfiletype.ext}",
desc=f"Writing geometry view file {i + 1}/{len(gvs)}, "
f"{gv.filename.name}{gv.vtkfiletype.ext}",
ncols=get_terminal_width() - 1,
file=sys.stdout,
disable=not config.sim_config.general["progressbars"],
@@ -96,7 +98,9 @@ class GeometryView:
def set_filename(self):
"""Constructs filename from user-supplied name and model run number."""
parts = config.get_model_config().output_file_path.parts
self.filename = Path(*parts[:-1], self.filenamebase + config.get_model_config().appendmodelnumber)
self.filename = Path(
*parts[:-1], self.filenamebase + config.get_model_config().appendmodelnumber
)
class GeometryViewLines(GeometryView):
@@ -122,7 +126,12 @@ class GeometryViewLines(GeometryView):
):
# Require contiguous for evtk library
ID = np.ascontiguousarray(
self.grid.ID[:, self.xs : self.xf : self.dx, self.ys : self.yf : self.dy, self.zs : self.zf : self.dz]
self.grid.ID[
:,
self.xs : self.xf : self.dx,
self.ys : self.yf : self.dy,
self.zs : self.zf : self.dz,
]
)
else:
# This array is contiguous by design
@@ -158,7 +167,15 @@ class GeometryViewLines(GeometryView):
offsets_size = np.arange(start=2, step=2, stop=len(x) + 1, dtype="int32").nbytes
connect_size = len(x) * np.dtype("int32").itemsize
cell_type_size = len(x) * np.dtype("uint8").itemsize
self.nbytes = x.nbytes + y.nbytes + z.nbytes + lines.nbytes + offsets_size + connect_size + cell_type_size
self.nbytes = (
x.nbytes
+ y.nbytes
+ z.nbytes
+ lines.nbytes
+ offsets_size
+ connect_size
+ cell_type_size
)
vtk_data = {"x": x, "y": y, "z": z, "data": lines, "comments": comments}
@@ -205,7 +222,11 @@ class GeometryViewVoxels(GeometryView):
):
# Require contiguous for evtk library
solid = np.ascontiguousarray(
self.grid.solid[self.xs : self.xf : self.dx, self.ys : self.yf : self.dy, self.zs : self.zf : self.dz]
self.grid.solid[
self.xs : self.xf : self.dx,
self.ys : self.yf : self.dy,
self.zs : self.zf : self.dz,
]
)
else:
# This array is contiguous by design
@@ -366,12 +387,18 @@ class GeometryObjects:
self.filename_materials = self.filename_materials.with_suffix(".txt")
# Sizes of arrays to write necessary to update progress bar
self.solidsize = (self.nx + 1) * (self.ny + 1) * (self.nz + 1) * np.dtype(np.uint32).itemsize
self.rigidsize = 18 * (self.nx + 1) * (self.ny + 1) * (self.nz + 1) * np.dtype(np.int8).itemsize
self.IDsize = 6 * (self.nx + 1) * (self.ny + 1) * (self.nz + 1) * np.dtype(np.uint32).itemsize
self.solidsize = (
(self.nx + 1) * (self.ny + 1) * (self.nz + 1) * np.dtype(np.uint32).itemsize
)
self.rigidsize = (
18 * (self.nx + 1) * (self.ny + 1) * (self.nz + 1) * np.dtype(np.int8).itemsize
)
self.IDsize = (
6 * (self.nx + 1) * (self.ny + 1) * (self.nz + 1) * np.dtype(np.uint32).itemsize
)
self.datawritesize = self.solidsize + self.rigidsize + self.IDsize
def write_hdf5(self, G, pbar):
def write_hdf5(self, title: str, G: FDTDGrid, pbar: tqdm):
"""Writes a geometry objects file in HDF5 format.
Args:
@@ -381,20 +408,34 @@ class GeometryObjects:
with h5py.File(self.filename_hdf5, "w") as fdata:
fdata.attrs["gprMax"] = __version__
fdata.attrs["Title"] = G.title
fdata.attrs["Title"] = title
fdata.attrs["dx_dy_dz"] = (G.dx, G.dy, G.dz)
# Get minimum and maximum integers of materials in geometry objects volume
minmat = np.amin(G.ID[:, self.xs : self.xf + 1, self.ys : self.yf + 1, self.zs : self.zf + 1])
maxmat = np.amax(G.ID[:, self.xs : self.xf + 1, self.ys : self.yf + 1, self.zs : self.zf + 1])
minmat = np.amin(
G.ID[:, self.xs : self.xf + 1, self.ys : self.yf + 1, self.zs : self.zf + 1]
)
maxmat = np.amax(
G.ID[:, self.xs : self.xf + 1, self.ys : self.yf + 1, self.zs : self.zf + 1]
)
fdata["/data"] = (
G.solid[self.xs : self.xf + 1, self.ys : self.yf + 1, self.zs : self.zf + 1].astype("int16") - minmat
G.solid[self.xs : self.xf + 1, self.ys : self.yf + 1, self.zs : self.zf + 1].astype(
"int16"
)
- minmat
)
pbar.update(self.solidsize)
fdata["/rigidE"] = G.rigidE[:, self.xs : self.xf + 1, self.ys : self.yf + 1, self.zs : self.zf + 1]
fdata["/rigidH"] = G.rigidH[:, self.xs : self.xf + 1, self.ys : self.yf + 1, self.zs : self.zf + 1]
fdata["/rigidE"] = G.rigidE[
:, self.xs : self.xf + 1, self.ys : self.yf + 1, self.zs : self.zf + 1
]
fdata["/rigidH"] = G.rigidH[
:, self.xs : self.xf + 1, self.ys : self.yf + 1, self.zs : self.zf + 1
]
pbar.update(self.rigidsize)
fdata["/ID"] = G.ID[:, self.xs : self.xf + 1, self.ys : self.yf + 1, self.zs : self.zf + 1] - minmat
fdata["/ID"] = (
G.ID[:, self.xs : self.xf + 1, self.ys : self.yf + 1, self.zs : self.zf + 1]
- minmat
)
pbar.update(self.IDsize)
# Write materials list to a text file
@@ -411,7 +452,9 @@ class GeometryObjects:
if "debye" in material.type:
dispersionstr = "#add_dispersion_debye: " f"{material.poles:g} "
for pole in range(material.poles):
dispersionstr += f"{material.deltaer[pole]:g} " f"{material.tau[pole]:g} "
dispersionstr += (
f"{material.deltaer[pole]:g} " f"{material.tau[pole]:g} "
)
elif "lorenz" in material.type:
dispersionstr = f"#add_dispersion_lorenz: " f"{material.poles:g} "
for pole in range(material.poles):
@@ -423,6 +466,8 @@ class GeometryObjects:
elif "drude" in material.type:
dispersionstr = f"#add_dispersion_drude: " f"{material.poles:g} "
for pole in range(material.poles):
dispersionstr += f"{material.tau[pole]:g} " f"{material.alpha[pole]:g} "
dispersionstr += (
f"{material.tau[pole]:g} " f"{material.alpha[pole]:g} "
)
dispersionstr += material.ID
fmaterials.write(dispersionstr + "\n")

580
gprMax/grid/fdtd_grid.py
查看文件

@@ -21,25 +21,22 @@ import itertools
import logging
import sys
from collections import OrderedDict
from typing import Iterable, List, Union
from typing import Any, Iterable, List, Union
import humanize
import numpy as np
from terminaltables import SingleTable
from tqdm import tqdm
from gprMax import config
from gprMax.cython.yee_cell_build import build_electric_components, build_magnetic_components
# from gprMax.geometry_outputs import GeometryObjects, GeometryView
from gprMax.materials import Material, process_materials
from gprMax.fractals import FractalVolume
from gprMax.materials import ListMaterial, Material, PeplinskiSoil, RangeMaterial, process_materials
from gprMax.pml import CFS, PML, build_pml, print_pml_info
from gprMax.receivers import Rx
from gprMax.sources import HertzianDipole, MagneticDipole, Source, VoltageSource
# from gprMax.subgrids.grid import SubGridBaseGrid
from gprMax.utilities.host_info import mem_check_build_all, mem_check_run_all
from gprMax.snapshots import Snapshot
from gprMax.sources import HertzianDipole, MagneticDipole, Source, TransmissionLine, VoltageSource
from gprMax.utilities.utilities import fft_power, get_terminal_width, round_value
from gprMax.waveforms import Waveform
logger = logging.getLogger(__name__)
@@ -49,22 +46,41 @@ class FDTDGrid:
accessing regularly used parameters.
"""
IDlookup = {"Ex": 0, "Ey": 1, "Ez": 2, "Hx": 3, "Hy": 4, "Hz": 5}
def __init__(self):
self.title = ""
self.name = "main_grid"
self.mem_use = 0
self.nx = 0
self.ny = 0
self.nz = 0
self.dx = 0
self.dy = 0
self.dz = 0
self.dt = 0
self.dt_mod = None # Time step stability factor
self.iteration = 0 # Current iteration number
self.iterations = 0 # Total number of iterations
self.timewindow = 0
self.dl = np.ones(3, dtype=float)
self.dt = 0.0
# Field Arrays
self.Ex: np.ndarray[Any, np.dtype[np.single]]
self.Ey: np.ndarray[Any, np.dtype[np.single]]
self.Ez: np.ndarray[Any, np.dtype[np.single]]
self.Hx: np.ndarray[Any, np.dtype[np.single]]
self.Hy: np.ndarray[Any, np.dtype[np.single]]
self.Hz: np.ndarray[Any, np.dtype[np.single]]
# Dispersive Arrays
self.Tx: np.ndarray[Any, np.dtype[np.single]]
self.Ty: np.ndarray[Any, np.dtype[np.single]]
self.Tz: np.ndarray[Any, np.dtype[np.single]]
# Geometry Arrays
self.solid: np.ndarray[Any, np.dtype[np.uint32]]
self.rigidE: np.ndarray[Any, np.dtype[np.int8]]
self.rigidH: np.ndarray[Any, np.dtype[np.int8]]
self.ID: np.ndarray[Any, np.dtype[np.uint32]]
# Update Coefficient Arrays
self.updatecoeffsE: np.ndarray
self.updatecoeffsH: np.ndarray
self.updatecoeffsdispersive: np.ndarray
# PML parameters - set some defaults to use if not user provided
self.pmls = {}
@@ -78,126 +94,63 @@ class FDTDGrid:
# corrections will be different.
self.pmls["thickness"] = OrderedDict((key, 10) for key in PML.boundaryIDs)
# TODO: Add type information.
# Currently importing GeometryObjects, GeometryView, and
# SubGridBaseGrid cause cyclic dependencies
# Materials used by this grid
self.materials: List[Material] = []
self.mixingmodels = []
self.averagevolumeobjects = True
self.fractalvolumes = []
self.geometryviews = []
self.geometryobjectswrite = []
self.waveforms = []
self.mixingmodels: List[Union[PeplinskiSoil, RangeMaterial, ListMaterial]] = []
self.fractalvolumes: List[FractalVolume] = []
# Sources and receivers contained inside this grid
self.waveforms: List[Waveform] = []
self.voltagesources: List[VoltageSource] = []
self.hertziandipoles: List[HertzianDipole] = []
self.magneticdipoles: List[MagneticDipole] = []
self.transmissionlines = []
self.transmissionlines: List[TransmissionLine] = []
self.rxs: List[Rx] = []
self.srcsteps: List[float] = [0, 0, 0]
self.rxsteps: List[float] = [0, 0, 0]
self.snapshots = []
self.subgrids = []
self.snapshots: List[Snapshot] = []
self.averagevolumeobjects = True
@property
def dx(self) -> float:
return self.dl[0]
@dx.setter
def dx(self, value: float):
self.dl[0] = value
@property
def dy(self) -> float:
return self.dl[1]
@dy.setter
def dy(self, value: float):
self.dl[1] = value
@property
def dz(self) -> float:
return self.dl[2]
@dz.setter
def dz(self, value: float):
self.dl[2] = value
def build(self) -> None:
# Print info on any subgrids
for subgrid in self.subgrids:
subgrid.print_info()
# Combine available grids
grids = [self] + self.subgrids
# Check for dispersive materials (and specific type)
if config.get_model_config().materials["maxpoles"] != 0:
# TODO: This sets materials["drudelorentz"] based only the
# last grid/subgrid. Is that correct?
for grid in grids:
config.get_model_config().materials["drudelorentz"] = any(
[m for m in grid.materials if "drude" in m.type or "lorentz" in m.type]
)
# Set data type if any dispersive materials (must be done before memory checks)
config.get_model_config().set_dispersive_material_types()
# Check memory requirements to build model/scene (different to memory
# requirements to run model when FractalVolumes/FractalSurfaces are
# used as these can require significant additional memory)
total_mem_build, mem_strs_build = mem_check_build_all(grids)
# Check memory requirements to run model
total_mem_run, mem_strs_run = mem_check_run_all(grids)
if total_mem_build > total_mem_run:
logger.info(
f'\nMemory required (estimated): {" + ".join(mem_strs_build)} + '
f"~{humanize.naturalsize(config.get_model_config().mem_overhead)} "
f"overhead = {humanize.naturalsize(total_mem_build)}"
)
else:
logger.info(
f'\nMemory required (estimated): {" + ".join(mem_strs_run)} + '
f"~{humanize.naturalsize(config.get_model_config().mem_overhead)} "
f"overhead = {humanize.naturalsize(total_mem_run)}"
)
# Build grids
for grid in grids:
# Set default CFS parameter for PMLs if not user provided
if not grid.pmls["cfs"]:
grid.pmls["cfs"] = [CFS()]
logger.info(print_pml_info(grid))
if not all(value == 0 for value in grid.pmls["thickness"].values()):
grid._build_pmls()
if grid.averagevolumeobjects:
grid._build_components()
grid._tm_grid_update()
grid._update_voltage_source_materials()
grid.initialise_field_arrays()
grid.initialise_std_update_coeff_arrays()
if config.get_model_config().materials["maxpoles"] > 0:
grid.initialise_dispersive_arrays()
grid.initialise_dispersive_update_coeff_array()
grid._build_materials()
# Check to see if numerical dispersion might be a problem
results = dispersion_analysis(grid)
if results["error"]:
logger.warning(
f"\nNumerical dispersion analysis [{grid.name}] "
f"not carried out as {results['error']}"
)
elif results["N"] < config.get_model_config().numdispersion["mingridsampling"]:
logger.exception(
f"\nNon-physical wave propagation in [{grid.name}] "
f"detected. Material '{results['material'].ID}' "
f"has wavelength sampled by {results['N']} cells, "
f"less than required minimum for physical wave "
f"propagation. Maximum significant frequency "
f"estimated as {results['maxfreq']:g}Hz"
)
raise ValueError
elif (
results["deltavp"]
and np.abs(results["deltavp"])
> config.get_model_config().numdispersion["maxnumericaldisp"]
):
logger.warning(
f"\n[{grid.name}] has potentially significant "
f"numerical dispersion. Estimated largest physical "
f"phase-velocity error is {results['deltavp']:.2f}% "
f"in material '{results['material'].ID}' whose "
f"wavelength sampled by {results['N']} cells. "
f"Maximum significant frequency estimated as "
f"{results['maxfreq']:g}Hz"
)
elif results["deltavp"]:
logger.info(
f"\nNumerical dispersion analysis [{grid.name}]: "
f"estimated largest physical phase-velocity error is "
f"{results['deltavp']:.2f}% in material '{results['material'].ID}' "
f"whose wavelength sampled by {results['N']} cells. "
f"Maximum significant frequency estimated as "
f"{results['maxfreq']:g}Hz"
)
# Set default CFS parameter for PMLs if not user provided
if not self.pmls["cfs"]:
self.pmls["cfs"] = [CFS()]
logger.info(print_pml_info(self))
if not all(value == 0 for value in self.pmls["thickness"].values()):
self._build_pmls()
if self.averagevolumeobjects:
self._build_components()
self._tm_grid_update()
self._create_voltage_source_materials()
self.initialise_field_arrays()
self.initialise_std_update_coeff_arrays()
if config.get_model_config().materials["maxpoles"] > 0:
self.initialise_dispersive_arrays()
self.initialise_dispersive_update_coeff_array()
self._build_materials()
def _build_pmls(self) -> None:
pbar = tqdm(
@@ -238,7 +191,7 @@ class FDTDGrid:
elif config.get_model_config().mode == "2D TMz":
self.tmz()
def _update_voltage_source_materials(self):
def _create_voltage_source_materials(self):
# Process any voltage sources (that have resistance) to create a new
# material at the source location
for voltagesource in self.voltagesources:
@@ -256,36 +209,36 @@ class FDTDGrid:
logger.info(materialstable.table)
def _update_positions(
self, items: Iterable[Union[Source, Rx]], step_size: List[float], step_number: int
self, items: Iterable[Union[Source, Rx]], step_size: List[int], step_number: int
) -> None:
if step_size[0] != 0 or step_size[1] != 0 or step_size[2] != 0:
for item in items:
if step_number == 0:
if (
item.xcoord + self.srcsteps[0] * config.sim_config.model_end < 0
or item.xcoord + self.srcsteps[0] * config.sim_config.model_end > self.nx
or item.ycoord + self.srcsteps[1] * config.sim_config.model_end < 0
or item.ycoord + self.srcsteps[1] * config.sim_config.model_end > self.ny
or item.zcoord + self.srcsteps[2] * config.sim_config.model_end < 0
or item.zcoord + self.srcsteps[2] * config.sim_config.model_end > self.nz
item.xcoord + step_size[0] * config.sim_config.model_end < 0
or item.xcoord + step_size[0] * config.sim_config.model_end > self.nx
or item.ycoord + step_size[1] * config.sim_config.model_end < 0
or item.ycoord + step_size[1] * config.sim_config.model_end > self.ny
or item.zcoord + step_size[2] * config.sim_config.model_end < 0
or item.zcoord + step_size[2] * config.sim_config.model_end > self.nz
):
raise ValueError
item.xcoord = item.xcoordorigin + step_number * step_size[0]
item.ycoord = item.ycoordorigin + step_number * step_size[1]
item.zcoord = item.zcoordorigin + step_number * step_size[2]
def update_simple_source_positions(self, step: int = 0) -> None:
def update_simple_source_positions(self, step_size: List[int], step: int = 0) -> None:
try:
self._update_positions(
itertools.chain(self.hertziandipoles, self.magneticdipoles), self.srcsteps, step
itertools.chain(self.hertziandipoles, self.magneticdipoles), step_size, step
)
except ValueError as e:
logger.exception("Source(s) will be stepped to a position outside the domain.")
raise ValueError from e
def update_receiver_positions(self, step: int = 0) -> None:
def update_receiver_positions(self, step_size: List[int], step: int = 0) -> None:
try:
self._update_positions(self.rxs, self.rxsteps, step)
self._update_positions(self.rxs, step_size, step)
except ValueError as e:
logger.exception("Receiver(s) will be stepped to a position outside the domain.")
raise ValueError from e
@@ -322,6 +275,9 @@ class FDTDGrid:
else:
return False
def get_waveform_by_id(self, waveform_id: str) -> Waveform:
return next(waveform for waveform in self.waveforms if waveform.ID == waveform_id)
def initialise_geometry_arrays(self):
"""Initialise an array for volumetric material IDs (solid);
boolean arrays for specifying whether materials can have dielectric
@@ -333,7 +289,6 @@ class FDTDGrid:
self.rigidE = np.zeros((12, self.nx, self.ny, self.nz), dtype=np.int8)
self.rigidH = np.zeros((6, self.nx, self.ny, self.nz), dtype=np.int8)
self.ID = np.ones((6, self.nx + 1, self.ny + 1, self.nz + 1), dtype=np.uint32)
self.IDlookup = {"Ex": 0, "Ey": 1, "Ez": 2, "Hx": 3, "Hy": 4, "Hz": 5}
def initialise_field_arrays(self):
"""Initialise arrays for the electric and magnetic field components."""
@@ -556,138 +511,237 @@ class FDTDGrid:
# binary representation of floating point number.
self.dt = round_value(self.dt, decimalplaces=decimal.getcontext().prec - 1)
def calculate_Ix(self, x: int, y: int, z: int) -> float:
"""Calculates the x-component of current at a grid position.
def dispersion_analysis(G):
"""Analysis of numerical dispersion (Taflove et al, 2005, p112) -
worse case of maximum frequency and minimum wavelength
Args:
x: x coordinate of position in grid
y: y coordinate of position in grid
z: z coordinate of position in grid
"""
Args:
G: FDTDGrid class describing a grid in a model.
if y == 0 or z == 0:
Ix = 0
else:
Ix = self.dy * (self.Hy[x, y, z - 1] - self.Hy[x, y, z]) + self.dz * (
self.Hz[x, y, z] - self.Hz[x, y - 1, z]
)
Returns:
results: dict of results from dispersion analysis.
"""
return Ix
# deltavp: physical phase velocity error (percentage)
# N: grid sampling density
# material: material with maximum permittivity
# maxfreq: maximum significant frequency
# error: error message
results = {"deltavp": None, "N": None, "material": None, "maxfreq": [], "error": ""}
def calculate_Iy(self, x: int, y: int, z: int) -> float:
"""Calculates the y-component of current at a grid position.
# Find maximum significant frequency
if G.waveforms:
for waveform in G.waveforms:
if waveform.type in ["sine", "contsine"]:
results["maxfreq"].append(4 * waveform.freq)
Args:
x: x coordinate of position in grid
y: y coordinate of position in grid
z: z coordinate of position in grid
"""
elif waveform.type == "impulse":
results["error"] = "impulse waveform used."
if x == 0 or z == 0:
Iy = 0
else:
Iy = self.dx * (self.Hx[x, y, z] - self.Hx[x, y, z - 1]) + self.dz * (
self.Hz[x - 1, y, z] - self.Hz[x, y, z]
)
elif waveform.type == "user":
results["error"] = "user waveform detected."
return Iy
else:
# Time to analyse waveform - 4*pulse_width as using entire
# time window can result in demanding FFT
waveform.calculate_coefficients()
iterations = round_value(4 * waveform.chi / G.dt)
iterations = min(iterations, G.iterations)
waveformvalues = np.zeros(G.iterations)
for iteration in range(G.iterations):
waveformvalues[iteration] = waveform.calculate_value(iteration * G.dt, G.dt)
def calculate_Iz(self, x: int, y: int, z: int) -> float:
"""Calculates the y-component of current at a grid position.
# Ensure source waveform is not being overly truncated before attempting any FFT
if np.abs(waveformvalues[-1]) < np.abs(np.amax(waveformvalues)) / 100:
# FFT
freqs, power = fft_power(waveformvalues, G.dt)
# Get frequency for max power
freqmaxpower = np.where(np.isclose(power, 0))[0][0]
Args:
x: x coordinate of position in grid
y: y coordinate of position in grid
z: z coordinate of position in grid
"""
# Set maximum frequency to a threshold drop from maximum power, ignoring DC value
try:
freqthres = (
np.where(
power[freqmaxpower:]
< -config.get_model_config().numdispersion["highestfreqthres"]
)[0][0]
+ freqmaxpower
)
results["maxfreq"].append(freqs[freqthres])
except ValueError:
results["error"] = (
"unable to calculate maximum power "
+ "from waveform, most likely due to "
+ "undersampling."
)
if x == 0 or y == 0:
Iz = 0
else:
Iz = self.dx * (self.Hx[x, y - 1, z] - self.Hx[x, y, z]) + self.dy * (
self.Hy[x, y, z] - self.Hy[x - 1, y, z]
)
# Ignore case where someone is using a waveform with zero amplitude, i.e. on a receiver
elif waveform.amp == 0:
pass
return Iz
# If waveform is truncated don't do any further analysis
else:
results["error"] = (
"waveform does not fit within specified "
+ "time window and is therefore being truncated."
)
else:
results["error"] = "no waveform detected."
if results["maxfreq"]:
results["maxfreq"] = max(results["maxfreq"])
# Find minimum wavelength (material with maximum permittivity)
maxer = 0
matmaxer = ""
for x in G.materials:
if x.se != float("inf"):
er = x.er
# If there are dispersive materials calculate the complex
# relative permittivity at maximum frequency and take the real part
if x.__class__.__name__ == "DispersiveMaterial":
er = x.calculate_er(results["maxfreq"])
er = er.real
if er > maxer:
maxer = er
matmaxer = x.ID
results["material"] = next(x for x in G.materials if x.ID == matmaxer)
# Minimum velocity
minvelocity = config.c / np.sqrt(maxer)
# Minimum wavelength
minwavelength = minvelocity / results["maxfreq"]
# Maximum spatial step
if "3D" in config.get_model_config().mode:
delta = max(G.dx, G.dy, G.dz)
elif "2D" in config.get_model_config().mode:
if G.nx == 1:
delta = max(G.dy, G.dz)
elif G.ny == 1:
delta = max(G.dx, G.dz)
elif G.nz == 1:
delta = max(G.dx, G.dy)
# Courant stability factor
S = (config.c * G.dt) / delta
# Grid sampling density
results["N"] = minwavelength / delta
# Check grid sampling will result in physical wave propagation
if (
int(np.floor(results["N"]))
>= config.get_model_config().numdispersion["mingridsampling"]
def dispersion_analysis(self, iterations: int):
# Check to see if numerical dispersion might be a problem
results = self._dispersion_analysis(iterations)
if results["error"]:
logger.warning(
f"\nNumerical dispersion analysis [{self.name}] "
f"not carried out as {results['error']}"
)
elif results["N"] < config.get_model_config().numdispersion["mingridsampling"]:
logger.exception(
f"\nNon-physical wave propagation in [{self.name}] "
f"detected. Material '{results['material'].ID}' "
f"has wavelength sampled by {results['N']} cells, "
f"less than required minimum for physical wave "
f"propagation. Maximum significant frequency "
f"estimated as {results['maxfreq']:g}Hz"
)
raise ValueError
elif (
results["deltavp"]
and np.abs(results["deltavp"])
> config.get_model_config().numdispersion["maxnumericaldisp"]
):
# Numerical phase velocity
vp = np.pi / (results["N"] * np.arcsin((1 / S) * np.sin((np.pi * S) / results["N"])))
logger.warning(
f"\n[{self.name}] has potentially significant "
f"numerical dispersion. Estimated largest physical "
f"phase-velocity error is {results['deltavp']:.2f}% "
f"in material '{results['material'].ID}' whose "
f"wavelength sampled by {results['N']} cells. "
f"Maximum significant frequency estimated as "
f"{results['maxfreq']:g}Hz"
)
elif results["deltavp"]:
logger.info(
f"\nNumerical dispersion analysis [{self.name}]: "
f"estimated largest physical phase-velocity error is "
f"{results['deltavp']:.2f}% in material '{results['material'].ID}' "
f"whose wavelength sampled by {results['N']} cells. "
f"Maximum significant frequency estimated as "
f"{results['maxfreq']:g}Hz"
)
# Physical phase velocity error (percentage)
results["deltavp"] = (((vp * config.c) - config.c) / config.c) * 100
def _dispersion_analysis(self, iterations: int):
"""Analysis of numerical dispersion (Taflove et al, 2005, p112) -
worse case of maximum frequency and minimum wavelength
# Store rounded down value of grid sampling density
results["N"] = int(np.floor(results["N"]))
Args:
G: FDTDGrid class describing a grid in a model.
return results
Returns:
results: dict of results from dispersion analysis.
"""
# deltavp: physical phase velocity error (percentage)
# N: grid sampling density
# material: material with maximum permittivity
# maxfreq: maximum significant frequency
# error: error message
results = {"deltavp": None, "N": None, "material": None, "maxfreq": [], "error": ""}
# Find maximum significant frequency
if self.waveforms:
for waveform in self.waveforms:
if waveform.type in ["sine", "contsine"]:
results["maxfreq"].append(4 * waveform.freq)
elif waveform.type == "impulse":
results["error"] = "impulse waveform used."
elif waveform.type == "user":
results["error"] = "user waveform detected."
else:
# Time to analyse waveform - 4*pulse_width as using entire
# time window can result in demanding FFT
waveform.calculate_coefficients()
iterations = round_value(4 * waveform.chi / self.dt)
iterations = min(iterations, iterations)
waveformvalues = np.zeros(iterations)
for iteration in range(iterations):
waveformvalues[iteration] = waveform.calculate_value(
iteration * self.dt, self.dt
)
# Ensure source waveform is not being overly truncated before attempting any FFT
if np.abs(waveformvalues[-1]) < np.abs(np.amax(waveformvalues)) / 100:
# FFT
freqs, power = fft_power(waveformvalues, self.dt)
# Get frequency for max power
freqmaxpower = np.where(np.isclose(power, 0))[0][0]
# Set maximum frequency to a threshold drop from maximum power, ignoring DC value
try:
freqthres = (
np.where(
power[freqmaxpower:]
< -config.get_model_config().numdispersion["highestfreqthres"]
)[0][0]
+ freqmaxpower
)
results["maxfreq"].append(freqs[freqthres])
except ValueError:
results["error"] = (
"unable to calculate maximum power "
+ "from waveform, most likely due to "
+ "undersampling."
)
# Ignore case where someone is using a waveform with zero amplitude, i.e. on a receiver
elif waveform.amp == 0:
pass
# If waveform is truncated don't do any further analysis
else:
results["error"] = (
"waveform does not fit within specified "
+ "time window and is therefore being truncated."
)
else:
results["error"] = "no waveform detected."
if results["maxfreq"]:
results["maxfreq"] = max(results["maxfreq"])
# Find minimum wavelength (material with maximum permittivity)
maxer = 0
matmaxer = ""
for x in self.materials:
if x.se != float("inf"):
er = x.er
# If there are dispersive materials calculate the complex
# relative permittivity at maximum frequency and take the real part
if x.__class__.__name__ == "DispersiveMaterial":
er = x.calculate_er(results["maxfreq"])
er = er.real
if er > maxer:
maxer = er
matmaxer = x.ID
results["material"] = next(x for x in self.materials if x.ID == matmaxer)
# Minimum velocity
minvelocity = config.c / np.sqrt(maxer)
# Minimum wavelength
minwavelength = minvelocity / results["maxfreq"]
# Maximum spatial step
if "3D" in config.get_model_config().mode:
delta = max(self.dx, self.dy, self.dz)
elif "2D" in config.get_model_config().mode:
if self.nx == 1:
delta = max(self.dy, self.dz)
elif self.ny == 1:
delta = max(self.dx, self.dz)
elif self.nz == 1:
delta = max(self.dx, self.dy)
# Courant stability factor
S = (config.c * self.dt) / delta
# Grid sampling density
results["N"] = minwavelength / delta
# Check grid sampling will result in physical wave propagation
if (
int(np.floor(results["N"]))
>= config.get_model_config().numdispersion["mingridsampling"]
):
# Numerical phase velocity
vp = np.pi / (
results["N"] * np.arcsin((1 / S) * np.sin((np.pi * S) / results["N"]))
)
# Physical phase velocity error (percentage)
results["deltavp"] = (((vp * config.c) - config.c) / config.c) * 100
# Store rounded down value of grid sampling density
results["N"] = int(np.floor(results["N"]))
return results

187
gprMax/model.py
查看文件

@@ -17,9 +17,9 @@
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
import datetime
import itertools
import logging
import sys
from typing import List
import humanize
import numpy as np
@@ -28,20 +28,17 @@ from colorama import Fore, Style, init
from gprMax.grid.cuda_grid import CUDAGrid
from gprMax.grid.opencl_grid import OpenCLGrid
from gprMax.subgrids.grid import SubGridBaseGrid
init()
from terminaltables import SingleTable
from tqdm import tqdm
import gprMax.config as config
from .cython.yee_cell_build import build_electric_components, build_magnetic_components
from .fields_outputs import write_hdf5_outputfile
from .geometry_outputs import save_geometry_views
from .grid.fdtd_grid import FDTDGrid, dispersion_analysis
from .materials import process_materials
from .pml import CFS, build_pml, print_pml_info
from .geometry_outputs import GeometryObjects, GeometryView, save_geometry_views
from .grid.fdtd_grid import FDTDGrid
from .snapshots import save_snapshots
from .utilities.host_info import mem_check_build_all, mem_check_run_all, set_omp_threads
from .utilities.utilities import get_terminal_width
@@ -53,7 +50,27 @@ class Model:
"""Builds and runs (solves) a model."""
def __init__(self):
self.title = ""
self.gnx = 0
self.gny = 0
self.gnz = 0
self.dt_mod = 1.0 # Time step stability factor
self.iteration = 0 # Current iteration number
self.iterations = 0 # Total number of iterations
self.timewindow = 0.0
self.srcsteps: List[int] = [0, 0, 0]
self.rxsteps: List[int] = [0, 0, 0]
self.G = self._create_grid()
self.subgrids: List[SubGridBaseGrid] = []
self.geometryviews: List[GeometryView] = []
self.geometryobjects: List[GeometryObjects] = []
# Monitor memory usage
self.p = None
@@ -63,6 +80,70 @@ class Model:
# later for use with CPU solver.
config.get_model_config().ompthreads = set_omp_threads(config.get_model_config().ompthreads)
@property
def nx(self) -> int:
return self.G.nx
@nx.setter
def nx(self, value: int):
self.G.nx = value
@property
def ny(self) -> int:
return self.G.ny
@ny.setter
def ny(self, value: int):
self.G.ny = value
@property
def nz(self) -> int:
return self.G.nz
@nz.setter
def nz(self, value: int):
self.G.nz = value
@property
def dx(self) -> float:
return self.G.dl[0]
@dx.setter
def dx(self, value: float):
self.G.dl[0] = value
@property
def dy(self) -> float:
return self.G.dl[0]
@dy.setter
def dy(self, value: float):
self.G.dl[1] = value
@property
def dz(self) -> float:
return self.G.dl[0]
@dz.setter
def dz(self, value: float):
self.G.dl[2] = value
@property
def dl(self) -> np.ndarray:
return self.G.dl
@dl.setter
def dl(self, value: np.ndarray):
self.G.dl = value
@property
def dt(self) -> float:
return self.G.dt
@dt.setter
def dt(self, value: float):
self.G.dt = value
def _create_grid(self) -> FDTDGrid:
"""Create grid object according to solver.
@@ -95,36 +176,92 @@ class Model:
# Adjust position of simple sources and receivers if required
model_num = config.sim_config.current_model
G.update_simple_source_positions(step=model_num)
G.update_receiver_positions(step=model_num)
G.update_simple_source_positions(self.srcsteps, step=model_num)
G.update_receiver_positions(self.rxsteps, step=model_num)
# Write files for any geometry views and geometry object outputs
gvs = G.geometryviews + [gv for sg in G.subgrids for gv in sg.geometryviews]
if not gvs and not G.geometryobjectswrite and config.sim_config.args.geometry_only:
if (
not self.geometryviews
and not self.geometryobjects
and config.sim_config.args.geometry_only
):
logger.exception("\nNo geometry views or geometry objects found.")
raise ValueError
save_geometry_views(gvs)
save_geometry_views(self.geometryviews)
if G.geometryobjectswrite:
if self.geometryobjects:
logger.info("")
for i, go in enumerate(G.geometryobjectswrite):
for i, go in enumerate(self.geometryobjects):
pbar = tqdm(
total=go.datawritesize,
unit="byte",
unit_scale=True,
desc=f"Writing geometry object file {i + 1}/{len(G.geometryobjectswrite)}, "
desc=f"Writing geometry object file {i + 1}/{len(self.geometryobjects)}, "
+ f"{go.filename_hdf5.name}",
ncols=get_terminal_width() - 1,
file=sys.stdout,
disable=not config.sim_config.general["progressbars"],
)
go.write_hdf5(G, pbar)
go.write_hdf5(self.title, self.G, pbar)
pbar.close()
logger.info("")
def build_geometry(self):
logger.info(config.get_model_config().inputfilestr)
self.G.build()
# Print info on any subgrids
for subgrid in self.subgrids:
subgrid.print_info()
# Combine available grids
grids = [self.G] + self.subgrids
self._check_for_dispersive_materials(grids)
self._check_memory_requirements(grids)
# TODO: Make this correctly set local nx, ny and nz when using MPI (likely use a function inside FDTDGrid/MPIGrid)
self.G.nx = self.gnx
self.G.ny = self.gny
self.G.nz = self.gnz
for grid in grids:
grid.build()
grid.dispersion_analysis(self.iterations)
def _check_for_dispersive_materials(self, grids: List[FDTDGrid]):
# Check for dispersive materials (and specific type)
if config.get_model_config().materials["maxpoles"] != 0:
# TODO: This sets materials["drudelorentz"] based only the
# last grid/subgrid. Is that correct?
for grid in grids:
config.get_model_config().materials["drudelorentz"] = any(
[m for m in grid.materials if "drude" in m.type or "lorentz" in m.type]
)
# Set data type if any dispersive materials (must be done before memory checks)
config.get_model_config().set_dispersive_material_types()
def _check_memory_requirements(self, grids: List[FDTDGrid]):
# Check memory requirements to build model/scene (different to memory
# requirements to run model when FractalVolumes/FractalSurfaces are
# used as these can require significant additional memory)
total_mem_build, mem_strs_build = mem_check_build_all(grids)
# Check memory requirements to run model
total_mem_run, mem_strs_run = mem_check_run_all(grids)
if total_mem_build > total_mem_run:
logger.info(
f'\nMemory required (estimated): {" + ".join(mem_strs_build)} + '
f"~{humanize.naturalsize(config.get_model_config().mem_overhead)} "
f"overhead = {humanize.naturalsize(total_mem_build)}"
)
else:
logger.info(
f'\nMemory required (estimated): {" + ".join(mem_strs_run)} + '
f"~{humanize.naturalsize(config.get_model_config().mem_overhead)} "
f"overhead = {humanize.naturalsize(total_mem_run)}"
)
def reuse_geometry(self):
s = (
@@ -136,8 +273,8 @@ class Model:
Fore.GREEN + f"{s} {'-' * (get_terminal_width() - 1 - len(s))}\n" + Style.RESET_ALL
)
logger.basic(config.get_model_config().inputfilestr)
for grid in [self.G] + self.G.subgrids:
grid.iteration = 0 # Reset current iteration number
self.iteration = 0 # Reset current iteration number
for grid in [self.G] + self.subgrids:
grid.reset_fields()
def write_output_data(self):
@@ -146,13 +283,13 @@ class Model:
"""
# Write output data to file if they are any receivers in any grids
sg_rxs = [True for sg in self.G.subgrids if sg.rxs]
sg_tls = [True for sg in self.G.subgrids if sg.transmissionlines]
sg_rxs = [True for sg in self.subgrids if sg.rxs]
sg_tls = [True for sg in self.subgrids if sg.transmissionlines]
if self.G.rxs or sg_rxs or self.G.transmissionlines or sg_tls:
write_hdf5_outputfile(config.get_model_config().output_file_path_ext, self.G)
write_hdf5_outputfile(config.get_model_config().output_file_path_ext, self.title, self)
# Write any snapshots to file for each grid
for grid in [self.G] + self.G.subgrids:
for grid in [self.G] + self.subgrids:
if grid.snapshots:
save_snapshots(grid)
@@ -204,14 +341,14 @@ class Model:
# Prepare iterator
if config.sim_config.general["progressbars"]:
iterator = tqdm(
range(self.G.iterations),
range(self.iterations),
desc="|--->",
ncols=get_terminal_width() - 1,
file=sys.stdout,
disable=not config.sim_config.general["progressbars"],
)
else:
iterator = range(self.G.iterations)
iterator = range(self.iterations)
# Run solver
solver.solve(iterator)

17
gprMax/receivers.py
查看文件

@@ -29,14 +29,14 @@ class Rx:
allowableoutputs_dev = allowableoutputs[:-3]
def __init__(self):
self.ID = None
self.ID: str
self.outputs = {}
self.xcoord = None
self.ycoord = None
self.zcoord = None
self.xcoordorigin = None
self.ycoordorigin = None
self.zcoordorigin = None
self.xcoord: int
self.ycoord: int
self.zcoord: int
self.xcoordorigin: int
self.ycoordorigin: int
self.zcoordorigin: int
def htod_rx_arrays(G, queue=None):
@@ -63,7 +63,8 @@ def htod_rx_arrays(G, queue=None):
# Array to store field components for receivers on compute device -
# rows are field components; columns are iterations; pages are receivers
rxs = np.zeros(
(len(Rx.allowableoutputs_dev), G.iterations, len(G.rxs)), dtype=config.sim_config.dtypes["float_or_double"]
(len(Rx.allowableoutputs_dev), G.iterations, len(G.rxs)),
dtype=config.sim_config.dtypes["float_or_double"],
)
# Copy arrays to compute device

91
gprMax/scene.py
查看文件

@@ -15,10 +15,10 @@
#
# You should have received a copy of the GNU General Public License
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
import logging
from typing import List, Union
from typing import List, Optional, Union
from gprMax import config
from gprMax.cmds_geometry.add_grass import AddGrass
from gprMax.cmds_geometry.add_surface_roughness import AddSurfaceRoughness
from gprMax.cmds_geometry.add_surface_water import AddSurfaceWater
@@ -28,9 +28,10 @@ from gprMax.cmds_multiuse import UserObjectMulti
from gprMax.cmds_singleuse import Discretisation, Domain, TimeWindow, UserObjectSingle
from gprMax.grid.fdtd_grid import FDTDGrid
from gprMax.materials import create_built_in_materials
from gprMax.model import Model
from gprMax.subgrids.grid import SubGridBaseGrid
from gprMax.subgrids.user_objects import SubGridBase as SubGridUserBase
from gprMax.user_inputs import create_user_input_points
from gprMax.user_inputs import MainGridUserInput, SubgridUserInput
logger = logging.getLogger(__name__)
@@ -39,6 +40,9 @@ class Scene:
"""Scene stores all of the user created objects."""
def __init__(self):
self.multiple_cmds: List[UserObjectMulti] = []
self.single_cmds: List[UserObjectSingle] = []
self.geometry_cmds: List[UserObjectGeometry] = []
self.multiple_cmds: List[UserObjectMulti] = []
self.single_cmds: List[UserObjectSingle] = []
self.geometry_cmds: List[UserObjectGeometry] = []
@@ -61,8 +65,8 @@ class Scene:
logger.exception("This object is unknown to gprMax")
raise ValueError
def build_obj(self, obj, grid):
"""Builds objects.
def build_grid_obj(self, obj: UserObjectGeometry, grid: FDTDGrid):
"""Builds objects in FDTDGrids.
Args:
obj: user object
@@ -75,13 +79,38 @@ class Scene:
logger.exception("Error creating user input object")
raise
def process_subgrid_cmds(self):
def build_model_obj(
self,
obj: Union[UserObjectSingle, UserObjectMulti],
model: Model,
subgrid: Optional[FDTDGrid] = None,
):
"""Builds objects in models.
Args:
obj: user object
model: Model being built
"""
grid = model.G if subgrid is None else subgrid
uip = create_user_input_points(grid, obj)
try:
obj.build(model, uip)
except ValueError:
logger.exception("Error creating user input object")
raise
def process_subgrid_cmds(self, model: Model):
"""Process all commands in any sub-grids."""
# Subgrid user objects
subgrid_cmds = [
sg_cmd for sg_cmd in self.multiple_cmds if isinstance(sg_cmd, SubGridUserBase)
]
subgrid_cmds = [
sg_cmd for sg_cmd in self.multiple_cmds if isinstance(sg_cmd, SubGridUserBase)
]
# Iterate through the user command objects under the subgrid user object
for sg_cmd in subgrid_cmds:
@@ -89,16 +118,19 @@ class Scene:
# object. This reference allows the multi and geo user objects
# to build in the correct subgrid.
sg = sg_cmd.subgrid
self.process_cmds(sg_cmd.children_multiple, sg)
self.process_cmds(sg_cmd.children_multiple, model, sg)
self.process_geocmds(sg_cmd.children_geometry, sg)
def process_cmds(
self, commands: Union[List[UserObjectMulti], List[UserObjectSingle]], grid: FDTDGrid
self,
commands: Union[List[UserObjectMulti], List[UserObjectSingle]],
model: Model,
subgrid: Optional[SubGridBaseGrid] = None,
):
"""Process list of commands."""
cmds_sorted = sorted(commands, key=lambda cmd: cmd.order)
for obj in cmds_sorted:
self.build_obj(obj, grid)
self.build_model_obj(obj, model, subgrid)
return self
@@ -107,7 +139,7 @@ class Scene:
proc_cmds = []
for obj in commands:
if isinstance(obj, (FractalBox, AddGrass, AddSurfaceRoughness, AddSurfaceWater)):
self.build_obj(obj, grid)
self.build_grid_obj(obj, grid)
if isinstance(obj, (FractalBox)):
proc_cmds.append(obj)
else:
@@ -115,11 +147,11 @@ class Scene:
# Process all geometry commands
for obj in proc_cmds:
self.build_obj(obj, grid)
self.build_grid_obj(obj, grid)
return self
def process_singlecmds(self, G: FDTDGrid):
def process_singlecmds(self, model: Model):
# Check for duplicate commands and warn user if they exist
cmds_unique = list(set(self.single_cmds))
if len(cmds_unique) != len(self.single_cmds):
@@ -137,29 +169,48 @@ class Scene:
)
raise ValueError
self.process_cmds(cmds_unique, G)
self.process_cmds(cmds_unique, model)
def create_internal_objects(self, G):
def create_internal_objects(self, model: Model):
"""Calls the UserObject.build() function in the correct way - API
presents the user with UserObjects in order to build the internal
Rx(), Cylinder() etc... objects.
"""
# Create pre-defined (built-in) materials
create_built_in_materials(G)
create_built_in_materials(model.G)
# Process commands that can only have a single instance
self.process_singlecmds(G)
self.process_singlecmds(model)
# Process main grid multiple commands
self.process_cmds(self.multiple_cmds, G)
self.process_cmds(self.multiple_cmds, model)
# Initialise geometry arrays for main and subgrids
for grid in [G] + G.subgrids:
for grid in [model.G] + model.subgrids:
grid.initialise_geometry_arrays()
# Process the main grid geometry commands
self.process_geocmds(self.geometry_cmds, G)
self.process_geocmds(self.geometry_cmds, model.G)
# Process all the commands for subgrids
self.process_subgrid_cmds()
self.process_subgrid_cmds(model)
def create_user_input_points(
grid: FDTDGrid, user_obj: Union[UserObjectSingle, UserObjectMulti, UserObjectGeometry]
) -> Union[MainGridUserInput, SubgridUserInput]:
"""Returns a point checker class based on the grid supplied."""
if isinstance(grid, SubGridBaseGrid):
# Local object configuration trumps. User can turn off autotranslate for
# specific objects.
if not user_obj.autotranslate and config.sim_config.args.autotranslate:
return MainGridUserInput(grid)
if config.sim_config.args.autotranslate:
return SubgridUserInput(grid)
else:
return MainGridUserInput(grid)
else:
return MainGridUserInput(grid)

90
gprMax/snapshots.py
查看文件

@@ -19,6 +19,7 @@
import logging
import sys
from pathlib import Path
from typing import Dict
import h5py
import numpy as np
@@ -55,7 +56,9 @@ def save_snapshots(grid):
leave=True,
unit="byte",
unit_scale=True,
desc=f"Writing snapshot file {i + 1} " f"of {len(grid.snapshots)}, " f"{snap.filename.name}",
desc=f"Writing snapshot file {i + 1} "
f"of {len(grid.snapshots)}, "
f"{snap.filename.name}",
ncols=get_terminal_width() - 1,
file=sys.stdout,
disable=not config.sim_config.general["progressbars"],
@@ -86,19 +89,19 @@ class Snapshot:
def __init__(
self,
xs=None,
ys=None,
zs=None,
xf=None,
yf=None,
zf=None,
dx=None,
dy=None,
dz=None,
time=None,
filename=None,
fileext=None,
outputs=None,
xs: int,
ys: int,
zs: int,
xf: int,
yf: int,
zf: int,
dx: int,
dy: int,
dz: int,
time: int,
filename: str,
fileext: str,
outputs: Dict[str, bool],
):
"""
Args:
@@ -142,7 +145,9 @@ class Snapshot:
else:
# If output is not required for snapshot just use a mimimal
# size of array - still required to pass to Cython function
self.snapfields[k] = np.zeros((1, 1, 1), dtype=config.sim_config.dtypes["float_or_double"])
self.snapfields[k] = np.zeros(
(1, 1, 1), dtype=config.sim_config.dtypes["float_or_double"]
)
def store(self, G):
"""Store (in memory) electric and magnetic field values for snapshot.
@@ -186,7 +191,7 @@ class Snapshot:
self.snapfields["Hz"],
)
def write_file(self, pbar, G):
def write_file(self, pbar: tqdm, G):
"""Writes snapshot file either as VTK ImageData (.vti) format
or HDF5 format (.h5) files
@@ -200,7 +205,7 @@ class Snapshot:
elif self.fileext == ".h5":
self.write_hdf5(pbar, G)
def write_vtk(self, pbar, G):
def write_vtk(self, pbar: tqdm, G):
"""Writes snapshot file in VTK ImageData (.vti) format.
Args:
@@ -208,11 +213,19 @@ class Snapshot:
G: FDTDGrid class describing a grid in a model.
"""
celldata = {k: self.snapfields[k] for k in ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"] if self.outputs.get(k)}
celldata = {
k: self.snapfields[k]
for k in ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"]
if self.outputs.get(k)
}
imageToVTK(
str(self.filename.with_suffix("")),
origin=((self.xs * self.dx * G.dx), (self.ys * self.dy * G.dy), (self.zs * self.dz * G.dz)),
origin=(
(self.xs * self.dx * G.dx),
(self.ys * self.dy * G.dy),
(self.zs * self.dz * G.dz),
),
spacing=((self.dx * G.dx), (self.dy * G.dy), (self.dz * G.dz)),
cellData=celldata,
)
@@ -225,7 +238,7 @@ class Snapshot:
* np.dtype(config.sim_config.dtypes["float_or_double"]).itemsize
)
def write_hdf5(self, pbar, G):
def write_hdf5(self, pbar: tqdm, G):
"""Writes snapshot file in HDF5 (.h5) format.
Args:
@@ -235,7 +248,8 @@ class Snapshot:
f = h5py.File(self.filename, "w")
f.attrs["gprMax"] = __version__
f.attrs["Title"] = G.title
# TODO: Output model name (title) and grid name? in snapshot output
# f.attrs["Title"] = G.title
f.attrs["nx_ny_nz"] = (self.nx, self.ny, self.nz)
f.attrs["dx_dy_dz"] = (self.dx * G.dx, self.dy * G.dy, self.dz * G.dz)
f.attrs["time"] = self.time * G.dt
@@ -271,35 +285,49 @@ def htod_snapshot_array(G, queue=None):
if config.sim_config.general["solver"] == "cuda":
# Blocks per grid - according to largest requested snapshot
Snapshot.bpg = (
int(np.ceil(((Snapshot.nx_max) * (Snapshot.ny_max) * (Snapshot.nz_max)) / Snapshot.tpb[0])),
int(
np.ceil(
((Snapshot.nx_max) * (Snapshot.ny_max) * (Snapshot.nz_max)) / Snapshot.tpb[0]
)
),
1,
1,
)
elif config.sim_config.general["solver"] == "opencl":
# Workgroup size - according to largest requested snapshot
Snapshot.wgs = (int(np.ceil(((Snapshot.nx_max) * (Snapshot.ny_max) * (Snapshot.nz_max)))), 1, 1)
Snapshot.wgs = (
int(np.ceil(((Snapshot.nx_max) * (Snapshot.ny_max) * (Snapshot.nz_max)))),
1,
1,
)
# 4D arrays to store snapshots on GPU, e.g. snapEx(time, x, y, z);
# if snapshots are not being stored on the GPU during the simulation then
# they are copied back to the host after each iteration, hence numsnaps = 1
numsnaps = 1 if config.get_model_config().device["snapsgpu2cpu"] else len(G.snapshots)
snapEx = np.zeros(
(numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=config.sim_config.dtypes["float_or_double"]
(numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
dtype=config.sim_config.dtypes["float_or_double"],
)
snapEy = np.zeros(
(numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=config.sim_config.dtypes["float_or_double"]
(numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
dtype=config.sim_config.dtypes["float_or_double"],
)
snapEz = np.zeros(
(numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=config.sim_config.dtypes["float_or_double"]
(numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
dtype=config.sim_config.dtypes["float_or_double"],
)
snapHx = np.zeros(
(numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=config.sim_config.dtypes["float_or_double"]
(numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
dtype=config.sim_config.dtypes["float_or_double"],
)
snapHy = np.zeros(
(numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=config.sim_config.dtypes["float_or_double"]
(numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
dtype=config.sim_config.dtypes["float_or_double"],
)
snapHz = np.zeros(
(numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=config.sim_config.dtypes["float_or_double"]
(numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
dtype=config.sim_config.dtypes["float_or_double"],
)
# Copy arrays to compute device
@@ -326,7 +354,9 @@ def htod_snapshot_array(G, queue=None):
return snapEx_dev, snapEy_dev, snapEz_dev, snapHx_dev, snapHy_dev, snapHz_dev
def dtoh_snapshot_array(snapEx_dev, snapEy_dev, snapEz_dev, snapHx_dev, snapHy_dev, snapHz_dev, i, snap):
def dtoh_snapshot_array(
snapEx_dev, snapEy_dev, snapEz_dev, snapHx_dev, snapHy_dev, snapHz_dev, i, snap
):
"""Copies snapshot array used on compute device back to snapshot objects and
store in format for Paraview.

16
gprMax/solvers.py
查看文件

@@ -17,6 +17,7 @@
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
import gprMax.config as config
from gprMax.model import Model
from .grid.cuda_grid import CUDAGrid
from .grid.fdtd_grid import FDTDGrid
@@ -53,16 +54,17 @@ class Solver:
self.updates.time_start()
for iteration in iterator:
self.updates.store_outputs()
self.updates.store_outputs(iteration)
self.updates.store_snapshots(iteration)
self.updates.update_magnetic()
self.updates.update_magnetic_pml()
self.updates.update_magnetic_sources()
self.updates.update_magnetic_sources(iteration)
if isinstance(self.updates, SubgridUpdates):
self.updates.hsg_2()
self.updates.update_electric_a()
self.updates.update_electric_pml()
self.updates.update_electric_sources()
self.updates.update_electric_sources(iteration)
# TODO: Increment iteration here if add Model to Solver
if isinstance(self.updates, SubgridUpdates):
self.updates.hsg_1()
self.updates.update_electric_b()
@@ -74,7 +76,7 @@ class Solver:
self.updates.cleanup()
def create_solver(grid: FDTDGrid) -> Solver:
def create_solver(model: Model) -> Solver:
"""Create configured solver object.
N.B. A large range of different functions exist to advance the time
@@ -86,14 +88,14 @@ def create_solver(grid: FDTDGrid) -> Solver:
substitution at runtime.
Args:
G: FDTDGrid class describing a grid in a model.
model: model containing the main grid and subgrids.
Returns:
solver: Solver object.
"""
grid = model.G
if config.sim_config.general["subgrid"]:
updates = create_subgrid_updates(grid)
updates = create_subgrid_updates(model)
if config.get_model_config().materials["maxpoles"] != 0:
# Set dispersive update functions for both SubgridUpdates and
# SubgridUpdaters subclasses

89
gprMax/sources.py
查看文件

@@ -21,8 +21,8 @@ from copy import deepcopy
import numpy as np
import gprMax.config as config
from gprMax.waveforms import Waveform
from .fields_outputs import Ix, Iy, Iz
from .utilities.utilities import round_value
@@ -30,40 +30,44 @@ class Source:
"""Super-class which describes a generic source."""
def __init__(self):
self.ID = None
self.ID: str
self.polarisation = None
self.xcoord = None
self.ycoord = None
self.zcoord = None
self.xcoordorigin = None
self.ycoordorigin = None
self.zcoordorigin = None
self.start = None
self.stop = None
self.waveformID = None
self.xcoord: int
self.ycoord: int
self.zcoord: int
self.xcoordorigin: int
self.ycoordorigin: int
self.zcoordorigin: int
self.start = 0.0
self.stop = 0.0
self.waveform: Waveform
def calculate_waveform_values(self, G):
def calculate_waveform_values(self, iterations: int, dt: float):
"""Calculates all waveform values for source for duration of simulation.
Args:
G: FDTDGrid class describing a grid in a model.
"""
# Waveform values for sources that need to be calculated on whole timesteps
self.waveformvalues_wholedt = np.zeros((G.iterations), dtype=config.sim_config.dtypes["float_or_double"])
self.waveformvalues_wholedt = np.zeros(
iterations, dtype=config.sim_config.dtypes["float_or_double"]
)
# Waveform values for sources that need to be calculated on half timesteps
self.waveformvalues_halfdt = np.zeros((G.iterations), dtype=config.sim_config.dtypes["float_or_double"])
self.waveformvalues_halfdt = np.zeros(
iterations, dtype=config.sim_config.dtypes["float_or_double"]
)
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
for iteration in range(G.iterations):
time = G.dt * iteration
for iteration in range(iterations):
time = dt * iteration
if time >= self.start and time <= self.stop:
# Set the time of the waveform evaluation to account for any
# delay in the start
time -= self.start
self.waveformvalues_wholedt[iteration] = waveform.calculate_value(time, G.dt)
self.waveformvalues_halfdt[iteration] = waveform.calculate_value(time + 0.5 * G.dt, G.dt)
self.waveformvalues_wholedt[iteration] = self.waveform.calculate_value(time, dt)
self.waveformvalues_halfdt[iteration] = self.waveform.calculate_value(
time + 0.5 * dt, dt
)
class VoltageSource(Source):
@@ -166,7 +170,7 @@ class HertzianDipole(Source):
def __init__(self):
super().__init__()
self.dl = None
self.dl = 0.0
def update_electric(self, iteration, updatecoeffsE, ID, Ex, Ey, Ez, G):
"""Updates electric field values for a Hertzian dipole.
@@ -274,7 +278,9 @@ def htod_src_arrays(sources, G, queue=None):
srcinfo1 = np.zeros((len(sources), 4), dtype=np.int32)
srcinfo2 = np.zeros((len(sources)), dtype=config.sim_config.dtypes["float_or_double"])
srcwaves = np.zeros((len(sources), G.iterations), dtype=config.sim_config.dtypes["float_or_double"])
srcwaves = np.zeros(
(len(sources), G.iterations), dtype=config.sim_config.dtypes["float_or_double"]
)
for i, src in enumerate(sources):
srcinfo1[i, 0] = src.xcoord
srcinfo1[i, 1] = src.ycoord
@@ -324,14 +330,16 @@ class TransmissionLine(Source):
which is attached virtually to a grid cell.
"""
def __init__(self, G):
def __init__(self, iterations: int, dt: float):
"""
Args:
G: FDTDGrid class describing a grid in a model.
iterations: number of iterations
dt: time step of the grid
"""
super().__init__()
self.resistance = None
self.iterations = iterations
# Coefficients for ABC termination of end of the transmission line
self.abcv0 = 0
@@ -339,11 +347,11 @@ class TransmissionLine(Source):
# Spatial step of transmission line (N.B if the magic time step is
# used it results in instabilities for certain impedances)
self.dl = np.sqrt(3) * config.c * G.dt
self.dl = np.sqrt(3) * config.c * dt
# Number of cells in the transmission line (initially a long line to
# calculate incident voltage and current); consider putting ABCs/PML at end
self.nl = round_value(0.667 * G.iterations)
self.nl = round_value(0.667 * self.iterations)
# Cell position of the one-way injector excitation in the transmission line
self.srcpos = 5
@@ -353,10 +361,10 @@ class TransmissionLine(Source):
self.voltage = np.zeros(self.nl, dtype=config.sim_config.dtypes["float_or_double"])
self.current = np.zeros(self.nl, dtype=config.sim_config.dtypes["float_or_double"])
self.Vinc = np.zeros(G.iterations, dtype=config.sim_config.dtypes["float_or_double"])
self.Iinc = np.zeros(G.iterations, dtype=config.sim_config.dtypes["float_or_double"])
self.Vtotal = np.zeros(G.iterations, dtype=config.sim_config.dtypes["float_or_double"])
self.Itotal = np.zeros(G.iterations, dtype=config.sim_config.dtypes["float_or_double"])
self.Vinc = np.zeros(self.iterations, dtype=config.sim_config.dtypes["float_or_double"])
self.Iinc = np.zeros(self.iterations, dtype=config.sim_config.dtypes["float_or_double"])
self.Vtotal = np.zeros(self.iterations, dtype=config.sim_config.dtypes["float_or_double"])
self.Itotal = np.zeros(self.iterations, dtype=config.sim_config.dtypes["float_or_double"])
def calculate_incident_V_I(self, G):
"""Calculates the incident voltage and current with a long length
@@ -367,7 +375,7 @@ class TransmissionLine(Source):
G: FDTDGrid class describing a grid in a model.
"""
for iteration in range(G.iterations):
for iteration in range(self.iterations):
self.Iinc[iteration] = self.current[self.antpos]
self.Vinc[iteration] = self.voltage[self.antpos]
self.update_current(iteration, G)
@@ -399,11 +407,15 @@ class TransmissionLine(Source):
# Update all the voltage values along the line
self.voltage[1 : self.nl] -= (
self.resistance * (config.c * G.dt / self.dl) * (self.current[1 : self.nl] - self.current[0 : self.nl - 1])
self.resistance
* (config.c * G.dt / self.dl)
* (self.current[1 : self.nl] - self.current[0 : self.nl - 1])
)
# Update the voltage at the position of the one-way injector excitation
self.voltage[self.srcpos] += (config.c * G.dt / self.dl) * self.waveformvalues_halfdt[iteration]
self.voltage[self.srcpos] += (config.c * G.dt / self.dl) * self.waveformvalues_halfdt[
iteration
]
# Update ABC before updating current
self.update_abc(G)
@@ -425,7 +437,9 @@ class TransmissionLine(Source):
# Update the current one cell before the position of the one-way injector excitation
self.current[self.srcpos - 1] += (
(1 / self.resistance) * (config.c * G.dt / self.dl) * self.waveformvalues_wholedt[iteration]
(1 / self.resistance)
* (config.c * G.dt / self.dl)
* self.waveformvalues_wholedt[iteration]
)
def update_electric(self, iteration, updatecoeffsE, ID, Ex, Ey, Ez, G):
@@ -458,6 +472,7 @@ class TransmissionLine(Source):
elif self.polarisation == "z":
Ez[i, j, k] = -self.voltage[self.antpos] / G.dz
# TODO: Add type information (if can avoid circular dependency)
def update_magnetic(self, iteration, updatecoeffsH, ID, Hx, Hy, Hz, G):
"""Updates current value in transmission line from magnetic field values
in the main grid.
@@ -478,12 +493,12 @@ class TransmissionLine(Source):
k = self.zcoord
if self.polarisation == "x":
self.current[self.antpos] = Ix(i, j, k, G.Hx, G.Hy, G.Hz, G)
self.current[self.antpos] = G.calculate_Ix(i, j, k)
elif self.polarisation == "y":
self.current[self.antpos] = Iy(i, j, k, G.Hx, G.Hy, G.Hz, G)
self.current[self.antpos] = G.calculate_Iy(i, j, k)
elif self.polarisation == "z":
self.current[self.antpos] = Iz(i, j, k, G.Hx, G.Hy, G.Hz, G)
self.current[self.antpos] = G.calculate_Iz(i, j, k)
self.update_current(iteration, G)

2
gprMax/subgrids/grid.py
查看文件

@@ -36,6 +36,8 @@ class SubGridBaseGrid(FDTDGrid, ABC):
# Name of the grid
self.name = kwargs["id"]
self.parent_grid: FDTDGrid
self.iterations = 0
self.filter = kwargs["filter"]

32
gprMax/subgrids/updates.py
查看文件

@@ -18,6 +18,10 @@
import logging
from gprMax.grid.fdtd_grid import FDTDGrid
from gprMax.model import Model
from gprMax.subgrids.grid import SubGridBaseGrid
from ..updates.cpu_updates import CPUUpdates
from .precursor_nodes import PrecursorNodes, PrecursorNodesFiltered
from .subgrid_hsg import SubGridHSG
@@ -25,24 +29,24 @@ from .subgrid_hsg import SubGridHSG
logger = logging.getLogger(__name__)
def create_updates(G):
def create_updates(model: Model):
"""Return the solver for the given subgrids."""
updaters = []
for sg in G.subgrids:
for sg in model.subgrids:
sg_type = type(sg)
if sg_type == SubGridHSG and sg.filter:
precursors = PrecursorNodesFiltered(G, sg)
precursors = PrecursorNodesFiltered(model.G, sg)
elif sg_type == SubGridHSG:
precursors = PrecursorNodes(G, sg)
precursors = PrecursorNodes(model.G, sg)
else:
logger.exception(f"{str(sg)} is not a subgrid type")
raise ValueError
sgu = SubgridUpdater(sg, precursors, G)
sgu = SubgridUpdater(sg, precursors, model.G)
updaters.append(sgu)
updates = SubgridUpdates(G, updaters)
updates = SubgridUpdates(model.G, updaters)
return updates
@@ -64,13 +68,13 @@ class SubgridUpdates(CPUUpdates):
sg_updater.hsg_2()
class SubgridUpdater(CPUUpdates):
class SubgridUpdater(CPUUpdates[SubGridBaseGrid]):
"""Handles 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):
def __init__(self, subgrid: SubGridBaseGrid, precursors: PrecursorNodes, G: FDTDGrid):
"""
Args:
subgrid: SubGrid3d instance to be updated.
@@ -82,7 +86,17 @@ class SubgridUpdater(CPUUpdates):
super().__init__(subgrid)
self.precursors = precursors
self.G = G
self.source_iteration = 0
self.iteration = 0
def store_outputs(self):
return super().store_outputs(self.iteration)
def update_electric_sources(self):
super().update_electric_sources(self.iteration)
self.iteration += 1
def update_magnetic_sources(self):
return super().update_magnetic_sources(self.iteration)
def hsg_1(self):
"""First half of the subgrid update. Takes the time step up to the main

30
gprMax/subgrids/user_objects.py
查看文件

@@ -43,6 +43,8 @@ class SubGridBase(UserObjectMulti):
super().__init__(**kwargs)
self.children_multiple: List[UserObjectMulti] = []
self.children_geometry: List[UserObjectGeometry] = []
self.children_multiple: List[UserObjectMulti] = []
self.children_geometry: List[UserObjectGeometry] = []
def add(self, node: Union[UserObjectMulti, UserObjectGeometry]):
"""Adds other user objects. Geometry and multi only."""
@@ -86,24 +88,24 @@ class SubGridBase(UserObjectMulti):
sg.ny = 2 * sg.n_boundary_cells_y + sg.nwy
sg.nz = 2 * sg.n_boundary_cells_z + sg.nwz
def set_iterations(self, sg: SubGridBaseGrid, main: FDTDGrid):
def set_iterations(self, sg: SubGridBaseGrid, model: Model):
"""Sets number of iterations that will take place in the subgrid."""
sg.iterations = main.iterations * sg.ratio
sg.iterations = model.iterations * sg.ratio
def setup(self, sg: SubGridBaseGrid, grid: FDTDGrid, uip: MainGridUserInput):
def setup(self, sg: SubGridBaseGrid, model: Model, uip: MainGridUserInput):
""" "Common setup to both all subgrid types."""
p1 = self.kwargs["p1"]
p2 = self.kwargs["p2"]
p1, p2 = uip.check_box_points(p1, p2, self.__str__())
self.set_discretisation(sg, grid)
self.set_discretisation(sg, model.G)
# Set temporal discretisation including any inherited time step
# stability factor from the main grid
sg.calculate_dt()
if grid.dt_mod:
sg.dt = sg.dt * grid.dt_mod
if model.dt_mod:
sg.dt = sg.dt * model.dt_mod
# Set the indices related to the subgrids main grid placement
self.set_main_grid_indices(sg, uip, p1, p2)
@@ -121,29 +123,27 @@ class SubGridBase(UserObjectMulti):
self.set_working_region_cells(sg)
self.set_total_cells(sg)
self.set_iterations(sg, grid)
self.set_iterations(sg, model)
self.set_name(sg)
# Copy a reference for the main grid to the sub grid
sg.parent_grid = grid
sg.timewindow = grid.timewindow
sg.parent_grid = model.G
# Copy a subgrid reference to self so that children.build(grid, uip)
# can access the correct grid.
self.subgrid = sg
# Copy over built in materials
sg.materials = [copy(m) for m in grid.materials if m.type == "builtin"]
sg.materials = [copy(m) for m in model.G.materials if m.type == "builtin"]
# Don't mix and match different subgrid types
for sg_made in grid.subgrids:
for sg_made in model.subgrids:
if type(sg) != type(sg_made):
logger.exception(f"{self.__str__()} please only use one type of subgrid")
raise ValueError
# Reference the subgrid under the main grid to which it belongs
grid.subgrids.append(sg)
model.subgrids.append(sg)
class SubGridHSG(SubGridBase):
@@ -200,7 +200,7 @@ class SubGridHSG(SubGridBase):
self.order = 18
self.hash = "#subgrid_hsg"
def build(self, grid: FDTDGrid, uip: MainGridUserInput) -> SubGridHSGUser:
def build(self, model: Model, uip: MainGridUserInput) -> SubGridHSGUser:
sg = SubGridHSGUser(**self.kwargs)
self.setup(sg, grid, uip)
self.setup(sg, model, uip)
return sg

22
gprMax/updates/cpu_updates.py
查看文件

@@ -17,6 +17,9 @@
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
from importlib import import_module
from typing import Generic
from typing_extensions import TypeVar
from gprMax import config
from gprMax.cython.fields_updates_normal import update_electric as update_electric_cpu
@@ -26,11 +29,13 @@ from gprMax.grid.fdtd_grid import FDTDGrid
from gprMax.updates.updates import Updates
from gprMax.utilities.utilities import timer
GridType = TypeVar("GridType", bound=FDTDGrid, default=FDTDGrid)
class CPUUpdates(Updates):
class CPUUpdates(Generic[GridType], Updates):
"""Defines update functions for CPU-based solver."""
def __init__(self, G: FDTDGrid):
def __init__(self, G: GridType):
"""
Args:
G: FDTDGrid class describing a grid in a model.
@@ -38,9 +43,9 @@ class CPUUpdates(Updates):
self.grid = G
def store_outputs(self):
def store_outputs(self, iteration):
"""Stores field component values for every receiver and transmission line."""
store_outputs_cpu(self.grid)
store_outputs_cpu(self.grid, iteration)
def store_snapshots(self, iteration):
"""Stores any snapshots.
@@ -74,11 +79,11 @@ class CPUUpdates(Updates):
for pml in self.grid.pmls["slabs"]:
pml.update_magnetic()
def update_magnetic_sources(self):
def update_magnetic_sources(self, iteration):
"""Updates magnetic field components from sources."""
for source in self.grid.transmissionlines + self.grid.magneticdipoles:
source.update_magnetic(
self.grid.iteration,
iteration,
self.grid.updatecoeffsH,
self.grid.ID,
self.grid.Hx,
@@ -134,7 +139,7 @@ class CPUUpdates(Updates):
for pml in self.grid.pmls["slabs"]:
pml.update_electric()
def update_electric_sources(self):
def update_electric_sources(self, iteration):
"""Updates electric field components from sources -
update any Hertzian dipole sources last.
"""
@@ -142,7 +147,7 @@ class CPUUpdates(Updates):
self.grid.voltagesources + self.grid.transmissionlines + self.grid.hertziandipoles
):
source.update_electric(
self.grid.iteration,
iteration,
self.grid.updatecoeffsE,
self.grid.ID,
self.grid.Ex,
@@ -150,7 +155,6 @@ class CPUUpdates(Updates):
self.grid.Ez,
self.grid,
)
self.grid.iteration += 1
def update_electric_b(self):
"""If there are any dispersive materials do 2nd part of dispersive

6
gprMax/updates/updates.py
查看文件

@@ -24,7 +24,7 @@ class Updates(ABC):
"""Defines update functions for a solver."""
@abstractmethod
def store_outputs(self) -> None:
def store_outputs(self, iteration: int) -> None:
"""Stores field component values for every receiver and transmission line."""
pass
@@ -48,7 +48,7 @@ class Updates(ABC):
pass
@abstractmethod
def update_magnetic_sources(self) -> None:
def update_magnetic_sources(self, iteration: int) -> None:
"""Updates magnetic field components from sources."""
pass
@@ -63,7 +63,7 @@ class Updates(ABC):
pass
@abstractmethod
def update_electric_sources(self) -> None:
def update_electric_sources(self, iteration: int) -> None:
"""Updates electric field components from sources -
update any Hertzian dipole sources last.
"""

48
gprMax/user_inputs.py
查看文件

@@ -15,14 +15,17 @@
#
# You should have received a copy of the GNU General Public License
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
from __future__ import annotations
import logging
from typing import Generic
import numpy as np
from typing_extensions import TypeVar
import gprMax.config as config
from gprMax.grid.fdtd_grid import FDTDGrid
from gprMax.subgrids.grid import SubGridBaseGrid
from .subgrids.grid import SubGridBaseGrid
from .utilities.utilities import round_value
logger = logging.getLogger(__name__)
@@ -36,28 +39,13 @@ logger = logging.getLogger(__name__)
encapulsated here.
"""
def create_user_input_points(grid, user_obj):
"""Returns a point checker class based on the grid supplied."""
if isinstance(grid, SubGridBaseGrid):
# Local object configuration trumps. User can turn off autotranslate for
# specific objects.
if not user_obj.autotranslate and config.sim_config.args.autotranslate:
return MainGridUserInput(grid)
if config.sim_config.args.autotranslate:
return SubgridUserInput(grid)
else:
return MainGridUserInput(grid)
else:
return MainGridUserInput(grid)
GridType = TypeVar("GridType", bound=FDTDGrid, default=FDTDGrid)
class UserInput:
class UserInput(Generic[GridType]):
"""Handles (x, y, z) points supplied by the user."""
def __init__(self, grid):
def __init__(self, grid: GridType):
self.grid = grid
def point_within_bounds(self, p, cmd_str, name):
@@ -92,7 +80,7 @@ class UserInput:
return p * self.grid.dl
class MainGridUserInput(UserInput):
class MainGridUserInput(UserInput[GridType]):
"""Handles (x, y, z) points supplied by the user in the main grid."""
def __init__(self, grid):
@@ -108,7 +96,9 @@ class MainGridUserInput(UserInput):
p = self.check_point(p, cmd_str, name)
if self.grid.within_pml(p):
logger.warning(f"'{cmd_str}' sources and receivers should not normally be positioned within the PML.")
logger.warning(
f"'{cmd_str}' sources and receivers should not normally be positioned within the PML."
)
return p
@@ -117,7 +107,9 @@ class MainGridUserInput(UserInput):
p2 = self.check_point(p2, cmd_str, name="upper")
if np.greater(p1, p2).any():
logger.exception(f"'{cmd_str}' the lower coordinates should be less than the upper coordinates.")
logger.exception(
f"'{cmd_str}' the lower coordinates should be less than the upper coordinates."
)
raise ValueError
return p1, p2
@@ -142,7 +134,7 @@ class MainGridUserInput(UserInput):
return super().discretise_point(p) * self.grid.dl
class SubgridUserInput(MainGridUserInput):
class SubgridUserInput(MainGridUserInput[SubGridBaseGrid]):
"""Handles (x, y, z) points supplied by the user in the subgrid.
This class autotranslates points from main grid to subgrid equivalent
(within IS). Useful if material traverse is not required.
@@ -152,7 +144,9 @@ class SubgridUserInput(MainGridUserInput):
super().__init__(grid)
# Defines the region exposed to the user
self.inner_bound = np.array([grid.n_boundary_cells_x, grid.n_boundary_cells_y, grid.n_boundary_cells_z])
self.inner_bound = np.array(
[grid.n_boundary_cells_x, grid.n_boundary_cells_y, grid.n_boundary_cells_z]
)
self.outer_bound = np.subtract([grid.nx, grid.ny, grid.nz], self.inner_bound)
@@ -186,7 +180,9 @@ class SubgridUserInput(MainGridUserInput):
# Provide user within a warning if they have placed objects within
# the OS non-working region.
if np.less(p_t, self.inner_bound).any() or np.greater(p_t, self.outer_bound).any():
logger.warning(f"'{cmd_str}' this object traverses the Outer Surface. This is an advanced feature.")
logger.warning(
f"'{cmd_str}' this object traverses the Outer Surface. This is an advanced feature."
)
return p_t
def discretise_static_point(self, p):

121
tests/grid/test_fdtd_grid.py 普通文件
查看文件

@@ -0,0 +1,121 @@
from typing import Callable
import numpy as np
import pytest
from numpy.testing import assert_allclose, assert_array_equal
from pytest import param
from gprMax.grid.fdtd_grid import FDTDGrid
def get_current_in_3d_grid(
get_current_func: Callable[[int, int, int], float], shape: tuple[int, ...]
) -> np.ndarray:
"""Helper function to get current as a 3D grid
Args:
get_current_func: Function that returns the current value at a
point on a grid.
shape: shape of the grid
Returns:
result: 3D grid containing current values
"""
result = np.empty(shape)
for i in range(shape[0]):
for j in range(shape[1]):
for k in range(shape[2]):
result[i, j, k] = get_current_func(i, j, k)
return result
def tile(a: float, b: float, n: int = 1) -> np.ndarray:
"""Helper function to tile two numbers
Args:
a: first value
b: second value
n: number of repititions
Returns:
c: tiled numpy array
"""
return np.tile([[[a, b], [b, a]], [[b, a], [a, b]]], (n, n, n))
@pytest.mark.parametrize(
"dy,dz,Hy,Hz,expected",
[
(0, 0, (1, 3), (0.5, 0.8), 0),
(0.1, 0.5, (0, 0), (0.5, 0.8), -0.15),
(0.1, 0.5, (1, 3), (0, 0), 0.2),
(0.1, 0, (1, 3), (0.5, 0.8), 0.2),
(0, 0.5, (1, 3), (0.5, 0.8), -0.15),
(0.1, 0.5, (1, 3), (0.5, 0.8), 0.05),
],
)
def test_calculate_Ix(dy, dz, Hy, Hz, expected, size=2):
grid = FDTDGrid()
grid.dy = dy
grid.dz = dz
grid.Hy = tile(Hy[0], Hy[1], size)
grid.Hz = tile(Hz[0], Hz[1], size)
actual_current = get_current_in_3d_grid(grid.calculate_Ix, grid.Hy.shape)
expected_current = tile(expected, -expected, size)
expected_current[:, 0, :] = 0
expected_current[:, :, 0] = 0
assert_allclose(actual_current, expected_current)
@pytest.mark.parametrize(
"dx,dz,Hx,Hz,expected",
[
(0, 0, (1, 3), (0.5, 0.8), 0),
(0.1, 0.5, (0, 0), (0.5, 0.8), -0.15),
(0.1, 0.5, (1, 3), (0, 0), 0.2),
(0.1, 0, (1, 3), (0.5, 0.8), 0.2),
(0, 0.5, (1, 3), (0.5, 0.8), -0.15),
(0.1, 0.5, (1, 3), (0.5, 0.8), 0.05),
],
)
def test_calculate_Iy(dx, dz, Hx, Hz, expected, size=2):
grid = FDTDGrid()
grid.dx = dx
grid.dz = dz
grid.Hx = tile(Hx[0], Hx[1], size)
grid.Hz = tile(Hz[0], Hz[1], size)
actual_current = get_current_in_3d_grid(grid.calculate_Iy, grid.Hx.shape)
expected_current = tile(-expected, expected, size)
expected_current[0, :, :] = 0
expected_current[:, :, 0] = 0
assert_allclose(actual_current, expected_current)
@pytest.mark.parametrize(
"dx,dy,Hx,Hy,expected",
[
(0, 0, (1, 3), (0.5, 0.8), 0),
(0.1, 0.5, (0, 0), (0.5, 0.8), -0.15),
(0.1, 0.5, (1, 3), (0, 0), 0.2),
(0.1, 0, (1, 3), (0.5, 0.8), 0.2),
(0, 0.5, (1, 3), (0.5, 0.8), -0.15),
(0.1, 0.5, (1, 3), (0.5, 0.8), 0.05),
],
)
def test_calculate_Iz(dx, dy, Hx, Hy, expected, size=2):
grid = FDTDGrid()
grid.dx = dx
grid.dy = dy
grid.Hx = tile(Hx[0], Hx[1], size)
grid.Hy = tile(Hy[0], Hy[1], size)
actual_current = get_current_in_3d_grid(grid.calculate_Iz, grid.Hx.shape)
expected_current = tile(expected, -expected, size)
expected_current[0, :, :] = 0
expected_current[:, 0, :] = 0
assert_allclose(actual_current, expected_current)