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

Merge branch '53-create-mpi-fractal-objects' into mpi

浏览代码
这个提交包含在：
Nathan Mannall
2025-05-22 15:58:25 +01:00
父节点 fa93dafcaf 0d0615c164
当前提交 6a754ca82f
共有 25 个文件被更改，包括 1556 次插入 和 641 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

29
gprMax/cython/fractals_generate.pyx
查看文件

@@ -26,6 +26,10 @@ from gprMax.config cimport float_or_double_complex
cpdef void generate_fractal2D(
int nx,
int ny,
int ox,
int oy,
int gx,
int gy,
int nthreads,
float D,
np.float64_t[:] weighting,
@@ -50,12 +54,16 @@ cpdef void generate_fractal2D(
cdef Py_ssize_t i, j
cdef double v2x, v2y, rr, B
cdef int sx, sy
sx = gx // 2
sy = gy // 2
for i in prange(nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(ny):
# Positional vector for current position
v2x = weighting[0] * i
v2y = weighting[1] * j
v2x = weighting[0] * ((i + ox + sx) % gx)
v2y = weighting[1] * ((j + oy + sy) % gy)
# Calulate norm of v2 - v1
rr = ((v2x - v1[0])**2 + (v2y - v1[1])**2)**(1/2)
@@ -70,6 +78,12 @@ cpdef void generate_fractal3D(
int nx,
int ny,
int nz,
int ox,
int oy,
int oz,
int gx,
int gy,
int gz,
int nthreads,
float D,
np.float64_t[:] weighting,
@@ -94,14 +108,19 @@ cpdef void generate_fractal3D(
cdef Py_ssize_t i, j, k
cdef double v2x, v2y, v2z, rr, B
cdef int sx, sy, sz
sx = gx // 2
sy = gy // 2
sz = gz // 2
for i in prange(nx, nogil=True, schedule='static', num_threads=nthreads):
for j in range(ny):
for k in range(nz):
# Positional vector for current position
v2x = weighting[0] * i
v2y = weighting[1] * j
v2z = weighting[2] * k
v2x = weighting[0] * ((i + ox + sx) % gx)
v2y = weighting[1] * ((j + oy + sy) % gy)
v2z = weighting[2] * ((k + oz + sz) % gz)
# Calulate norm of v2 - v1
rr = ((v2x - v1[0])**2 + (v2y - v1[1])**2 + (v2z - v1[2])**2)**(1/2)

340
gprMax/fractals.py
查看文件

@@ -1,340 +0,0 @@
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
#
# This file is part of gprMax.
#
# gprMax is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# gprMax is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
import numpy as np
from scipy import fftpack
import gprMax.config as config
from .cython.fractals_generate import generate_fractal2D, generate_fractal3D
from .utilities.utilities import round_value
np.seterr(divide="raise")
class FractalSurface:
"""Fractal surfaces."""
surfaceIDs = ["xminus", "xplus", "yminus", "yplus", "zminus", "zplus"]
def __init__(self, xs, xf, ys, yf, zs, zf, dimension, seed):
"""
Args:
xs, xf, ys, yf, zs, zf: floats for the extent of the fractal surface
(one pair of coordinates must be equal
to correctly define a surface).
dimension: float for the fractal dimension that controls the fractal
distribution.
seed: int for seed value for random number generator.
"""
self.ID = None
self.surfaceID = None
self.xs = xs
self.xf = xf
self.ys = ys
self.yf = yf
self.zs = zs
self.zf = zf
self.nx = xf - xs
self.ny = yf - ys
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.weighting = np.array([1, 1], dtype=np.float64)
self.fractalrange = (0, 0)
self.filldepth = 0
self.grass = []
def get_surface_dims(self):
"""Gets the dimensions of the fractal surface based on surface plane."""
if self.xs == self.xf:
surfacedims = (self.ny, self.nz)
elif self.ys == self.yf:
surfacedims = (self.nx, self.nz)
elif self.zs == self.zf:
surfacedims = (self.nx, self.ny)
return surfacedims
def generate_fractal_surface(self):
"""Generate a 2D array with a fractal distribution."""
surfacedims = self.get_surface_dims()
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,
]
)
# 2D array of random numbers to be convolved with the fractal function
rng = np.random.default_rng(seed=self.seed)
A = rng.standard_normal(size=(surfacedims[0], surfacedims[1]))
# 2D FFT
A = fftpack.fftn(A)
# Shift the zero frequency component to the centre of the array
A = fftpack.fftshift(A)
# Generate fractal
generate_fractal2D(
surfacedims[0],
surfacedims[1],
config.get_model_config().ompthreads,
self.dimension,
self.weighting,
v1,
A,
self.fractalsurface,
)
# Shift the zero frequency component to start of the array
self.fractalsurface = fftpack.ifftshift(self.fractalsurface)
# Set DC component of FFT to zero
self.fractalsurface[0, 0] = 0
# Take the real part (numerical errors can give rise to an imaginary part)
# of the IFFT, and convert type to floattype. N.B calculation of fractals
# must always be carried out at double precision, i.e. float64, complex128
self.fractalsurface = np.real(fftpack.ifftn(self.fractalsurface)).astype(
config.sim_config.dtypes["float_or_double"], copy=False
)
# Scale the fractal volume according to requested range
fractalmin = np.amin(self.fractalsurface)
fractalmax = np.amax(self.fractalsurface)
fractalrange = fractalmax - fractalmin
self.fractalsurface = (
self.fractalsurface * ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange)
+ self.fractalrange[0]
- ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * fractalmin
)
class FractalVolume:
"""Fractal volumes."""
def __init__(self, xs, xf, ys, yf, zs, zf, dimension, seed):
"""
Args:
xs, xf, ys, yf, zs, zf: floats for the extent of the fractal volume.
dimension: float for the fractal dimension that controls the fractal
distribution.
seed: int for seed value for random number generator.
"""
self.ID = None
self.operatingonID = None
self.xs = xs
self.xf = xf
self.ys = ys
self.yf = yf
self.zs = zs
self.zf = zf
self.nx = xf - xs
self.ny = yf - ys
self.nz = zf - zs
self.originalxs = xs
self.originalxf = xf
self.originalys = ys
self.originalyf = yf
self.originalzs = zs
self.originalzf = zf
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.weighting = np.array([1, 1, 1], dtype=np.float64)
self.nbins = 0
self.fractalsurfaces = []
def generate_fractal_volume(self):
"""Generate a 3D volume with a fractal distribution."""
# Scale filter according to size of fractal volume
if self.nx == 1:
filterscaling = np.amin(np.array([self.ny, self.nz])) / np.array([self.ny, self.nz])
filterscaling = np.insert(filterscaling, 0, 1)
elif self.ny == 1:
filterscaling = np.amin(np.array([self.nx, self.nz])) / np.array([self.nx, self.nz])
filterscaling = np.insert(filterscaling, 1, 1)
elif self.nz == 1:
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]
)
# Adjust weighting to account for filter scaling
self.weighting = np.multiply(self.weighting, filterscaling)
self.fractalvolume = np.zeros((self.nx, self.ny, self.nz), dtype=self.dtype)
# 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,
]
)
# 3D array of random numbers to be convolved with the fractal function
rng = np.random.default_rng(seed=self.seed)
A = rng.standard_normal(size=(self.nx, self.ny, self.nz))
# 3D FFT
A = fftpack.fftn(A)
# Shift the zero frequency component to the centre of the array
A = fftpack.fftshift(A)
# Generate fractal
generate_fractal3D(
self.nx,
self.ny,
self.nz,
config.get_model_config().ompthreads,
self.dimension,
self.weighting,
v1,
A,
self.fractalvolume,
)
# Shift the zero frequency component to the start of the array
self.fractalvolume = fftpack.ifftshift(self.fractalvolume)
# Set DC component of FFT to zero
self.fractalvolume[0, 0, 0] = 0
# Take the real part (numerical errors can give rise to an imaginary part)
# of the IFFT, and convert type to floattype. N.B calculation of fractals
# must always be carried out at double precision, i.e. float64, complex128
self.fractalvolume = np.real(fftpack.ifftn(self.fractalvolume)).astype(
config.sim_config.dtypes["float_or_double"], copy=False
)
# Bin fractal values
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
)
def generate_volume_mask(self):
"""Generate a 3D volume to use as a mask for adding rough surfaces,
water and grass/roots. Zero signifies the mask is not set, one
signifies the mask is set.
"""
self.mask = np.zeros((self.nx, self.ny, self.nz), dtype=np.int8)
maskxs = self.originalxs - self.xs
maskxf = (self.originalxf - self.originalxs) + maskxs
maskys = self.originalys - self.ys
maskyf = (self.originalyf - self.originalys) + maskys
maskzs = self.originalzs - self.zs
maskzf = (self.originalzf - self.originalzs) + maskzs
self.mask[maskxs:maskxf, maskys:maskyf, maskzs:maskzf] = 1
class Grass:
"""Geometry information for blades of grass."""
def __init__(self, numblades, seed):
"""
Args:
numblades: int for the number of blades of grass.
seed: int for seed value for random number generator.
"""
self.numblades = numblades
self.geometryparams = np.zeros(
(self.numblades, 6), dtype=config.sim_config.dtypes["float_or_double"]
)
self.seed = seed
self.set_geometry_parameters()
def set_geometry_parameters(self):
"""Sets randomly defined parameters that will be used to calculate
blade and root geometries.
"""
self.R1 = np.random.default_rng(seed=self.seed)
self.R2 = np.random.default_rng(seed=self.seed)
self.R3 = np.random.default_rng(seed=self.seed)
self.R4 = np.random.default_rng(seed=self.seed)
self.R5 = np.random.default_rng(seed=self.seed)
self.R6 = np.random.default_rng(seed=self.seed)
for i in range(self.numblades):
self.geometryparams[i, 0] = 10 + 20 * self.R1.random()
self.geometryparams[i, 1] = 10 + 20 * self.R2.random()
self.geometryparams[i, 2] = self.R3.choice([-1, 1])
self.geometryparams[i, 3] = self.R4.choice([-1, 1])
def calculate_blade_geometry(self, blade, height):
"""Calculates the x and y coordinates for a given height of grass blade.
Args:
blade: int for the numeric ID of grass blade.
height: float for the height of grass blade.
Returns:
x, y: floats for the x and y coordinates of grass blade.
"""
x = (
self.geometryparams[blade, 2]
* (height / self.geometryparams[blade, 0])
* (height / self.geometryparams[blade, 0])
)
y = (
self.geometryparams[blade, 3]
* (height / self.geometryparams[blade, 1])
* (height / self.geometryparams[blade, 1])
)
x = round_value(x)
y = round_value(y)
return x, y
def calculate_root_geometry(self, root, depth):
"""Calculates the x and y coordinates for a given depth of grass root.
Args:
root: int for the umeric ID of grass root.
depth: float for the depth of grass root.
Returns:
x, y: floats for the x and y coordinates of grass root.
"""
self.geometryparams[root, 4] += -1 + 2 * self.R5.random()
self.geometryparams[root, 5] += -1 + 2 * self.R6.random()
x = round(self.geometryparams[root, 4])
y = round(self.geometryparams[root, 5])
return x, y

428
gprMax/fractals/fractal_surface.py 普通文件
查看文件

@@ -0,0 +1,428 @@
import logging
from typing import List, Optional, Tuple
import numpy as np
import numpy.typing as npt
from mpi4py import MPI
from mpi4py_fft import PFFT, DistArray, newDistArray
from mpi4py_fft.pencil import Subcomm
from scipy import fftpack
from gprMax import config
from gprMax.cython.fractals_generate import generate_fractal2D
from gprMax.fractals.grass import Grass
from gprMax.fractals.mpi_utilities import calculate_starts_and_subshape, create_mpi_type
from gprMax.utilities.mpi import Dim, Dir, get_relative_neighbour
logger = logging.getLogger(__name__)
np.seterr(divide="raise")
class FractalSurface:
"""Fractal surfaces."""
surfaceIDs = ["xminus", "xplus", "yminus", "yplus", "zminus", "zplus"]
def __init__(self, xs, xf, ys, yf, zs, zf, dimension, seed):
"""
Args:
xs, xf, ys, yf, zs, zf: floats for the extent of the fractal surface
(one pair of coordinates must be equal
to correctly define a surface).
dimension: float for the fractal dimension that controls the fractal
distribution.
seed: int for seed value for random number generator.
"""
self.ID = None
self.surfaceID = None
self.start = np.array([xs, ys, zs], dtype=np.int32)
self.stop = np.array([xf, yf, zf], dtype=np.int32)
self.dtype = np.dtype(np.complex128)
self.seed = seed
self.dimension = (
dimension # Fractal dimension from: http://dx.doi.org/10.1017/CBO9781139174695
)
self.weighting = np.array([1, 1], dtype=np.float64)
self.fractalrange: Tuple[int, int] = (0, 0)
self.filldepth = 0
self.grass: List[Grass] = []
@property
def xs(self) -> int:
return self.start[0]
@xs.setter
def xs(self, value: int):
self.start[0] = value
@property
def ys(self) -> int:
return self.start[1]
@ys.setter
def ys(self, value: int):
self.start[1] = value
@property
def zs(self) -> int:
return self.start[2]
@zs.setter
def zs(self, value: int):
self.start[2] = value
@property
def xf(self) -> int:
return self.stop[0]
@xf.setter
def xf(self, value: int):
self.stop[0] = value
@property
def yf(self) -> int:
return self.stop[1]
@yf.setter
def yf(self, value: int):
self.stop[1] = value
@property
def zf(self) -> int:
return self.stop[2]
@zf.setter
def zf(self, value: int):
self.stop[2] = value
@property
def size(self) -> npt.NDArray[np.int32]:
return self.stop - self.start
@property
def nx(self) -> int:
return self.xf - self.xs
@property
def ny(self) -> int:
return self.yf - self.ys
@property
def nz(self) -> int:
return self.zf - self.zs
def get_surface_dims(self):
"""Gets the dimensions of the fractal surface based on surface plane."""
if self.xs == self.xf:
surfacedims = (self.ny, self.nz)
elif self.ys == self.yf:
surfacedims = (self.nx, self.nz)
elif self.zs == self.zf:
surfacedims = (self.nx, self.ny)
return surfacedims
def generate_fractal_surface(self):
"""Generate a 2D array with a fractal distribution."""
surfacedims = self.get_surface_dims()
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,
]
)
# 2D array of random numbers to be convolved with the fractal function
rng = np.random.default_rng(seed=self.seed)
A = rng.standard_normal(size=(surfacedims[0], surfacedims[1]))
# 2D FFT
A = fftpack.fftn(A)
# Generate fractal
generate_fractal2D(
surfacedims[0],
surfacedims[1],
0,
0,
surfacedims[0],
surfacedims[1],
config.get_model_config().ompthreads,
self.dimension,
self.weighting,
v1,
A,
self.fractalsurface,
)
# Set DC component of FFT to zero
self.fractalsurface[0, 0] = 0
# Take the real part (numerical errors can give rise to an imaginary part)
# of the IFFT, and convert type to floattype. N.B calculation of fractals
# must always be carried out at double precision, i.e. float64, complex128
self.fractalsurface = np.real(fftpack.ifftn(self.fractalsurface)).astype(
config.sim_config.dtypes["float_or_double"], copy=False
)
# Scale the fractal volume according to requested range
fractalmin = np.amin(self.fractalsurface)
fractalmax = np.amax(self.fractalsurface)
fractalrange = fractalmax - fractalmin
self.fractalsurface = (
self.fractalsurface * ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange)
+ self.fractalrange[0]
- ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * fractalmin
)
class MPIFractalSurface(FractalSurface):
def __init__(
self,
xs: int,
xf: int,
ys: int,
yf: int,
zs: int,
zf: int,
dimension: float,
seed: Optional[int],
comm: MPI.Cartcomm,
upper_bound: npt.NDArray[np.int32],
):
super().__init__(xs, xf, ys, yf, zs, zf, dimension, seed)
self.comm = comm
self.upper_bound = upper_bound
def generate_fractal_surface(self):
"""Generate a 2D array with a fractal distribution."""
if self.xs == self.xf:
color = self.xs
static_dimension = Dim.X
dims = [Dim.Y, Dim.Z]
elif self.ys == self.yf:
color = self.ys
static_dimension = Dim.Y
dims = [Dim.X, Dim.Z]
elif self.zs == self.zf:
color = self.zs
static_dimension = Dim.Z
dims = [Dim.X, Dim.Y]
# Exit early if this rank does not contain the Fractal Surface
if (
any(self.stop[dims] <= 0)
or any(self.start[dims] >= self.upper_bound[dims])
or self.fractalrange[1] <= 0
or self.fractalrange[0] >= self.upper_bound[static_dimension]
):
self.comm.Split(MPI.UNDEFINED)
# Update start and stop to local bounds
self.start = np.maximum(self.start, 0)
self.start = np.minimum(self.start, self.upper_bound)
self.stop = np.maximum(self.stop, 0)
self.stop = np.minimum(self.stop, self.upper_bound)
return
else:
# Create new cartsesian communicator for the Fractal Surface
comm = self.comm.Split(color=color)
assert isinstance(comm, MPI.Intracomm)
min_coord = np.array(self.comm.coords, dtype=np.int32)[dims]
max_coord = min_coord + 1
comm.Allreduce(MPI.IN_PLACE, min_coord, MPI.MIN)
comm.Allreduce(MPI.IN_PLACE, max_coord, MPI.MAX)
cart_dims = (max_coord - min_coord).tolist()
self.comm = comm.Create_cart(cart_dims)
# Check domain decomosition is valid for the Fractal Volume
if all([dim > 1 for dim in self.comm.dims]):
raise ValueError(
"Fractal surface must be positioned such that its MPI decomposition is 1 in at least"
f" 1 dimension. Current decompostion is: {self.comm.dims}"
)
surfacedims = self.get_surface_dims()
# Positional vector at centre of array, scaled by weighting
v1 = self.weighting * surfacedims / 2
subcomm = Subcomm(self.comm)
A = DistArray(self.size[dims], subcomm, dtype=self.dtype)
fft = PFFT(
None,
axes=tuple(np.argsort(self.comm.dims)[::-1]),
darray=A,
collapse=False,
backend="fftw",
)
# Decomposition of A may be different to the MPIGrid
A_shape = np.array(A.shape)
A_substart = np.array(A.substart)
# 3D array of random numbers to be convolved with the fractal function
rng = np.random.default_rng(seed=self.seed)
for index in np.ndindex(*A.global_shape):
index = np.array(index)
if any(index < A_substart) or any(index >= A_substart + A_shape):
rng.standard_normal()
else:
index -= A_substart
A[index[0], index[1]] = rng.standard_normal()
A_hat = newDistArray(fft)
assert isinstance(A_hat, DistArray)
# 2D FFT
fft.forward(A, A_hat, normalize=False)
# Generate fractal
generate_fractal2D(
A_hat.shape[0],
A_hat.shape[1],
A_hat.substart[0],
A_hat.substart[1],
A_hat.global_shape[0],
A_hat.global_shape[1],
config.get_model_config().ompthreads,
self.dimension,
self.weighting,
v1,
A_hat,
A_hat,
)
# Set DC component of FFT to zero
if all(A_substart == 0):
A_hat[0, 0] = 0
# Inverse 2D FFT transform
fft.backward(A_hat, A, normalize=True)
# Take the real part (numerical errors can give rise to an imaginary part)
# of the IFFT, and convert type to floattype. N.B calculation of fractals
# must always be carried out at double precision, i.e. float64, complex128
A = np.real(A).astype(config.sim_config.dtypes["float_or_double"], copy=False)
# Allreduce to get min and max values in the fractal surface
min_value = np.array(np.amin(A), dtype=config.sim_config.dtypes["float_or_double"])
max_value = np.array(np.amax(A), dtype=config.sim_config.dtypes["float_or_double"])
self.comm.Allreduce(MPI.IN_PLACE, min_value, MPI.MIN)
self.comm.Allreduce(MPI.IN_PLACE, max_value, MPI.MAX)
# Scale the fractal volume according to requested range
fractalrange = max_value - min_value
A = (
A * ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange)
+ self.fractalrange[0]
- ((self.fractalrange[1] - self.fractalrange[0]) / fractalrange) * min_value
)
# Distribute A (DistArray) to match the MPIGrid decomposition
local_shape = (np.minimum(self.stop, self.upper_bound) - np.maximum(self.start, 0))[dims]
self.fractalsurface = np.zeros(
local_shape,
dtype=config.sim_config.dtypes["float_or_double"],
)
# Negative means send to negative neighbour
# Positive means receive from negative neighbour
negative_offset = np.where(self.start[dims] >= 0, A_substart, self.start[dims] + A_substart)
# Negative means send to positive neighbour
# Positive means receive from positive neighbour
positive_offset = np.minimum(self.stop, self.upper_bound)[dims] - (
self.start[dims] + A_substart + A_shape
)
dirs = np.full(2, Dir.NONE)
starts, subshape = calculate_starts_and_subshape(
A_shape, -negative_offset, -positive_offset, dirs, sending=True
)
ends = starts + subshape
A_local = A[starts[0] : ends[0], starts[1] : ends[1]]
starts, subshape = calculate_starts_and_subshape(
local_shape, negative_offset, positive_offset, dirs
)
ends = starts + subshape
self.fractalsurface[starts[0] : ends[0], starts[1] : ends[1]] = A_local
requests: List[MPI.Request] = []
# Need to check neighbours in each direction (2D plane)
sections = [
(Dir.NEG, Dir.NONE),
(Dir.POS, Dir.NONE),
(Dir.NONE, Dir.NEG),
(Dir.NONE, Dir.POS),
(Dir.NEG, Dir.NEG),
(Dir.NEG, Dir.POS),
(Dir.POS, Dir.NEG),
(Dir.POS, Dir.POS),
]
for section in sections:
dirs[0] = section[0]
dirs[1] = section[1]
rank = get_relative_neighbour(self.comm, dirs)
# Skip if no neighbour
if rank == -1:
continue
# Check if any data to send
if all(
np.select(
[dirs == Dir.NEG, dirs == Dir.POS],
[negative_offset <= 0, positive_offset <= 0],
dirs == Dir.NONE,
)
):
mpi_type = create_mpi_type(
A_shape, -negative_offset, -positive_offset, dirs, sending=True
)
logger.debug(
f"Sending fractal surface to rank {rank}, MPI type={mpi_type.decode()}"
)
self.comm.Isend([A, mpi_type], rank)
# Check if any data to receive
if all(
np.select(
[dirs == Dir.NEG, dirs == Dir.POS],
[negative_offset > 0, positive_offset > 0],
dirs == Dir.NONE,
)
):
mpi_type = create_mpi_type(local_shape, negative_offset, positive_offset, dirs)
logger.debug(
f"Receiving fractal surface from rank {rank}, MPI type={mpi_type.decode()}"
)
request = self.comm.Irecv([self.fractalsurface, mpi_type], rank)
requests.append(request)
if len(requests) > 0:
requests[0].Waitall(requests)
# Update start and stop to local bounds
self.start = np.maximum(self.start, 0)
self.start = np.minimum(self.start, self.upper_bound)
self.stop = np.maximum(self.stop, 0)
self.stop = np.minimum(self.stop, self.upper_bound)
logger.debug(
f"Generated fractal surface: start={self.start}, stop={self.stop}, size={self.size}, fractalrange={self.fractalrange}"
)

524
gprMax/fractals/fractal_volume.py 普通文件
查看文件

@@ -0,0 +1,524 @@
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
#
# This file is part of gprMax.
#
# gprMax is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# gprMax is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with gprMax. If not, see <http://www.gnu.org/licenses/>.
import logging
from typing import List, Optional, Tuple, Union
import numpy as np
import numpy.typing as npt
from mpi4py import MPI
from mpi4py_fft import PFFT, DistArray, newDistArray
from mpi4py_fft.pencil import Subcomm
from scipy import fftpack
from gprMax import config
from gprMax.cython.fractals_generate import generate_fractal3D
from gprMax.fractals.fractal_surface import FractalSurface
from gprMax.fractals.mpi_utilities import calculate_starts_and_subshape, create_mpi_type
from gprMax.materials import ListMaterial, PeplinskiSoil, RangeMaterial
from gprMax.utilities.mpi import Dim, Dir, get_relative_neighbour
from gprMax.utilities.utilities import round_value
logger = logging.getLogger(__name__)
np.seterr(divide="raise")
class FractalVolume:
"""Fractal volumes."""
def __init__(
self,
xs: int,
xf: int,
ys: int,
yf: int,
zs: int,
zf: int,
dimension: float,
seed: Optional[int],
):
"""
Args:
xs, xf, ys, yf, zs, zf: floats for the extent of the fractal volume.
dimension: float for the fractal dimension that controls the fractal
distribution.
seed: int for seed value for random number generator.
"""
self.ID = None
self.operatingonID = None
self.start = np.array([xs, ys, zs], dtype=np.int32)
self.stop = np.array([xf, yf, zf], dtype=np.int32)
self.original_start = self.start.copy()
self.original_stop = self.stop.copy()
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.weighting = np.array([1, 1, 1], dtype=np.float64)
self.nbins = 0
self.mixingmodel: Optional[Union[PeplinskiSoil, RangeMaterial, ListMaterial]] = None
self.fractalsurfaces: List[FractalSurface] = []
@property
def xs(self) -> int:
return self.start[0]
@xs.setter
def xs(self, value: int):
self.start[0] = value
@property
def ys(self) -> int:
return self.start[1]
@ys.setter
def ys(self, value: int):
self.start[1] = value
@property
def zs(self) -> int:
return self.start[2]
@zs.setter
def zs(self, value: int):
self.start[2] = value
@property
def xf(self) -> int:
return self.stop[0]
@xf.setter
def xf(self, value: int):
self.stop[0] = value
@property
def yf(self) -> int:
return self.stop[1]
@yf.setter
def yf(self, value: int):
self.stop[1] = value
@property
def zf(self) -> int:
return self.stop[2]
@zf.setter
def zf(self, value: int):
self.stop[2] = value
@property
def size(self) -> npt.NDArray[np.int32]:
return self.stop - self.start
@property
def nx(self) -> int:
return self.xf - self.xs
@property
def ny(self) -> int:
return self.yf - self.ys
@property
def nz(self) -> int:
return self.zf - self.zs
@property
def originalxs(self) -> int:
return self.original_start[0]
@originalxs.setter
def originalxs(self, value: int):
self.original_start[0] = value
@property
def originalys(self) -> int:
return self.original_start[1]
@originalys.setter
def originalys(self, value: int):
self.original_start[1] = value
@property
def originalzs(self) -> int:
return self.original_start[2]
@originalzs.setter
def originalzs(self, value: int):
self.original_start[2] = value
@property
def originalxf(self) -> int:
return self.original_stop[0]
@originalxf.setter
def originalxf(self, value: int):
self.original_stop[0] = value
@property
def originalyf(self) -> int:
return self.original_stop[1]
@originalyf.setter
def originalyf(self, value: int):
self.original_stop[1] = value
@property
def originalzf(self) -> int:
return self.original_stop[2]
@originalzf.setter
def originalzf(self, value: int):
self.original_stop[2] = value
def generate_fractal_volume(self) -> bool:
"""Generate a 3D volume with a fractal distribution."""
# Scale filter according to size of fractal volume
if self.nx == 1:
filterscaling = np.amin(np.array([self.ny, self.nz])) / np.array([self.ny, self.nz])
filterscaling = np.insert(filterscaling, 0, 1)
elif self.ny == 1:
filterscaling = np.amin(np.array([self.nx, self.nz])) / np.array([self.nx, self.nz])
filterscaling = np.insert(filterscaling, 1, 1)
elif self.nz == 1:
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]
)
# Adjust weighting to account for filter scaling
self.weighting = np.multiply(self.weighting, filterscaling)
self.fractalvolume = np.zeros((self.nx, self.ny, self.nz), dtype=self.dtype)
# 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,
]
)
# 3D array of random numbers to be convolved with the fractal function
rng = np.random.default_rng(seed=self.seed)
A = rng.standard_normal(size=(self.nx, self.ny, self.nz))
# 3D FFT
A = fftpack.fftn(A)
# Generate fractal
generate_fractal3D(
self.nx,
self.ny,
self.nz,
0,
0,
0,
self.nx,
self.ny,
self.nz,
config.get_model_config().ompthreads,
self.dimension,
self.weighting,
v1,
A,
self.fractalvolume,
)
# Set DC component of FFT to zero
self.fractalvolume[0, 0, 0] = 0
# Take the real part (numerical errors can give rise to an imaginary part)
# of the IFFT, and convert type to floattype. N.B calculation of fractals
# must always be carried out at double precision, i.e. float64, complex128
self.fractalvolume = np.real(fftpack.ifftn(self.fractalvolume)).astype(
config.sim_config.dtypes["float_or_double"], copy=False
)
# Bin fractal values
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
)
return True
def generate_volume_mask(self):
"""Generate a 3D volume to use as a mask for adding rough surfaces,
water and grass/roots. Zero signifies the mask is not set, one
signifies the mask is set.
"""
self.mask = np.zeros((self.nx, self.ny, self.nz), dtype=np.int8)
maskxs = self.originalxs - self.xs
maskxf = (self.originalxf - self.originalxs) + maskxs
maskys = self.originalys - self.ys
maskyf = (self.originalyf - self.originalys) + maskys
maskzs = self.originalzs - self.zs
maskzf = (self.originalzf - self.originalzs) + maskzs
self.mask[maskxs:maskxf, maskys:maskyf, maskzs:maskzf] = 1
class MPIFractalVolume(FractalVolume):
def __init__(
self,
xs: int,
xf: int,
ys: int,
yf: int,
zs: int,
zf: int,
dimension: float,
seed: Optional[int],
comm: MPI.Cartcomm,
upper_bound: npt.NDArray[np.int32],
):
super().__init__(xs, xf, ys, yf, zs, zf, dimension, seed)
self.comm = comm
self.upper_bound = upper_bound
# Limit the original start and stop to within the local bounds
self.original_start = np.maximum(self.original_start, 0)
self.original_stop = np.minimum(self.original_stop, self.upper_bound)
# Ensure original_stop is not less than original_start
self.original_stop = np.where(
self.original_stop < self.original_start, self.original_start, self.original_stop
)
def generate_fractal_volume(self) -> bool:
"""Generate a 3D volume with a fractal distribution."""
# Exit early if this rank does not contain the Fractal Volume
# The size of a fractal volume can increase if a Fractal Surface
# is attached. Hence the check needs to happen here once that
# has happened.
if any(self.stop <= 0) or any(self.start >= self.upper_bound):
self.comm.Split(MPI.UNDEFINED)
return False
else:
# Create new cartsesian communicator for the Fractal Volume
comm = self.comm.Split()
assert isinstance(comm, MPI.Intracomm)
min_coord = np.array(self.comm.coords, dtype=np.int32)
max_coord = min_coord + 1
comm.Allreduce(MPI.IN_PLACE, min_coord, MPI.MIN)
comm.Allreduce(MPI.IN_PLACE, max_coord, MPI.MAX)
self.comm = comm.Create_cart((max_coord - min_coord).tolist())
# Check domain decomosition is valid for the Fractal Volume
if all([dim > 1 for dim in self.comm.dims]):
raise ValueError(
"Fractal volume must be positioned such that its MPI decomposition is 1 in at least"
f" 1 dimension. Current decompostion is: {self.comm.dims}"
)
# Scale filter according to size of fractal volume
sorted_size = np.sort(self.size)
min_size = sorted_size[1] if sorted_size[0] == 1 else sorted_size[0]
filterscaling = np.where(self.size == 1, 1, min_size / self.size)
# Adjust weighting to account for filter scaling
self.weighting = np.multiply(self.weighting, filterscaling)
# Positional vector at centre of array, scaled by weighting
v1 = self.weighting * self.size / 2
subcomm = Subcomm(self.comm)
A = DistArray(self.size, subcomm, dtype=self.dtype)
fft = PFFT(
None,
axes=tuple(np.argsort(self.comm.dims)[::-1]),
darray=A,
collapse=False,
backend="fftw",
)
# Decomposition of A may be different to the MPIGrid
A_shape = np.array(A.shape)
A_substart = np.array(A.substart)
static_dimension = Dim(A.alignment)
# 3D array of random numbers to be convolved with the fractal function
rng = np.random.default_rng(seed=self.seed)
for index in np.ndindex(*A.global_shape):
index = np.array(index)
if any(index < A_substart) or any(index >= A_substart + A_shape):
rng.standard_normal()
else:
index -= A_substart
A[index[0], index[1], index[2]] = rng.standard_normal()
A_hat = newDistArray(fft)
assert isinstance(A_hat, DistArray)
# 3D FFT
fft.forward(A, A_hat, normalize=False)
# Generate fractal
generate_fractal3D(
A_hat.shape[0],
A_hat.shape[1],
A_hat.shape[2],
A_hat.substart[0],
A_hat.substart[1],
A_hat.substart[2],
A_hat.global_shape[0],
A_hat.global_shape[1],
A_hat.global_shape[2],
config.get_model_config().ompthreads,
self.dimension,
self.weighting,
v1,
A_hat,
A_hat,
)
# Set DC component of FFT to zero
if all(A_substart == 0):
A_hat[0, 0, 0] = 0
# Inverse 3D FFT transform
fft.backward(A_hat, A, normalize=True)
# Take the real part (numerical errors can give rise to an imaginary part)
# of the IFFT, and convert type to floattype. N.B calculation of fractals
# must always be carried out at double precision, i.e. float64, complex128
A = np.real(A).astype(config.sim_config.dtypes["float_or_double"], copy=False)
# Allreduce to get min and max values in the fractal volume
min_value = np.array(np.amin(A), dtype=config.sim_config.dtypes["float_or_double"])
max_value = np.array(np.amax(A), dtype=config.sim_config.dtypes["float_or_double"])
self.comm.Allreduce(MPI.IN_PLACE, min_value, MPI.MIN)
self.comm.Allreduce(MPI.IN_PLACE, max_value, MPI.MAX)
# Bin fractal values
bins = np.linspace(min_value, max_value, self.nbins)
for j in range(A_shape[1]):
for k in range(A_shape[2]):
A[:, j, k] = np.digitize(A[:, j, k], bins, right=True)
# Distribute A (DistArray) to match the MPIGrid decomposition
local_shape = np.minimum(self.stop, self.upper_bound) - np.maximum(self.start, 0)
self.fractalvolume = np.zeros(
local_shape,
dtype=config.sim_config.dtypes["float_or_double"],
)
# Negative means send to negative neighbour
# Positive means receive from negative neighbour
negative_offset = np.where(self.start >= 0, A_substart, self.start + A_substart)
# Negative means send to positive neighbour
# Positive means receive from positive neighbour
positive_offset = np.minimum(self.stop, self.upper_bound) - (
self.start + A_substart + A_shape
)
dirs = np.full(3, Dir.NONE)
starts, subshape = calculate_starts_and_subshape(
A_shape, -negative_offset, -positive_offset, dirs, sending=True
)
ends = starts + subshape
A_local = A[starts[0] : ends[0], starts[1] : ends[1], starts[2] : ends[2]]
starts, subshape = calculate_starts_and_subshape(
local_shape, negative_offset, positive_offset, dirs
)
ends = starts + subshape
self.fractalvolume[starts[0] : ends[0], starts[1] : ends[1], starts[2] : ends[2]] = A_local
requests: List[MPI.Request] = []
# Need to check neighbours in each direction (2D plane)
sections = [
(Dir.NEG, Dir.NONE),
(Dir.POS, Dir.NONE),
(Dir.NONE, Dir.NEG),
(Dir.NONE, Dir.POS),
(Dir.NEG, Dir.NEG),
(Dir.NEG, Dir.POS),
(Dir.POS, Dir.NEG),
(Dir.POS, Dir.POS),
]
# Dimensions of the 2D plane
dims = [dim for dim in Dim if dim != static_dimension]
for section in sections:
dirs[dims[0]] = section[0]
dirs[dims[1]] = section[1]
rank = get_relative_neighbour(self.comm, dirs)
# Skip if no neighbour
if rank == -1:
continue
# Check if any data to send
if all(
np.select(
[dirs == Dir.NEG, dirs == Dir.POS],
[negative_offset <= 0, positive_offset <= 0],
dirs == Dir.NONE,
)
):
mpi_type = create_mpi_type(
A_shape, -negative_offset, -positive_offset, dirs, sending=True
)
logger.debug(f"Sending fractal volume to rank {rank}, MPI type={mpi_type.decode()}")
self.comm.Isend([A, mpi_type], rank)
# Check if any data to receive
if all(
np.select(
[dirs == Dir.NEG, dirs == Dir.POS],
[negative_offset > 0, positive_offset > 0],
dirs == Dir.NONE,
)
):
mpi_type = create_mpi_type(local_shape, negative_offset, positive_offset, dirs)
logger.debug(
f"Receiving fractal volume from rank {rank}, MPI type={mpi_type.decode()}"
)
request = self.comm.Irecv([self.fractalvolume, mpi_type], rank)
requests.append(request)
if len(requests) > 0:
requests[0].Waitall(requests)
# Update start and stop to local bounds
self.start = np.maximum(self.start, 0)
self.stop = np.minimum(self.stop, self.upper_bound)
logger.debug(
f"Generated fractal volume: start={self.start}, stop={self.stop}, size={self.size}"
)
return True

86
gprMax/fractals/grass.py 普通文件
查看文件

@@ -0,0 +1,86 @@
import numpy as np
from gprMax import config
from gprMax.utilities.utilities import round_value
np.seterr(divide="raise")
class Grass:
"""Geometry information for blades of grass."""
def __init__(self, numblades, seed):
"""
Args:
numblades: int for the number of blades of grass.
seed: int for seed value for random number generator.
"""
self.numblades = numblades
self.geometryparams = np.zeros(
(self.numblades, 6), dtype=config.sim_config.dtypes["float_or_double"]
)
self.seed = seed
self.set_geometry_parameters()
def set_geometry_parameters(self):
"""Sets randomly defined parameters that will be used to calculate
blade and root geometries.
"""
self.R1 = np.random.default_rng(seed=self.seed)
self.R2 = np.random.default_rng(seed=self.seed)
self.R3 = np.random.default_rng(seed=self.seed)
self.R4 = np.random.default_rng(seed=self.seed)
self.R5 = np.random.default_rng(seed=self.seed)
self.R6 = np.random.default_rng(seed=self.seed)
for i in range(self.numblades):
self.geometryparams[i, 0] = 10 + 20 * self.R1.random()
self.geometryparams[i, 1] = 10 + 20 * self.R2.random()
self.geometryparams[i, 2] = self.R3.choice([-1, 1])
self.geometryparams[i, 3] = self.R4.choice([-1, 1])
def calculate_blade_geometry(self, blade, height):
"""Calculates the x and y coordinates for a given height of grass blade.
Args:
blade: int for the numeric ID of grass blade.
height: float for the height of grass blade.
Returns:
x, y: floats for the x and y coordinates of grass blade.
"""
x = (
self.geometryparams[blade, 2]
* (height / self.geometryparams[blade, 0])
* (height / self.geometryparams[blade, 0])
)
y = (
self.geometryparams[blade, 3]
* (height / self.geometryparams[blade, 1])
* (height / self.geometryparams[blade, 1])
)
x = round_value(x)
y = round_value(y)
return x, y
def calculate_root_geometry(self, root, depth):
"""Calculates the x and y coordinates for a given depth of grass root.
Args:
root: int for the umeric ID of grass root.
depth: float for the depth of grass root.
Returns:
x, y: floats for the x and y coordinates of grass root.
"""
self.geometryparams[root, 4] += -1 + 2 * self.R5.random()
self.geometryparams[root, 5] += -1 + 2 * self.R6.random()
x = round(self.geometryparams[root, 4])
y = round(self.geometryparams[root, 5])
return x, y

62
gprMax/fractals/mpi_utilities.py 普通文件
查看文件

@@ -0,0 +1,62 @@
from typing import Tuple
import numpy as np
import numpy.typing as npt
from mpi4py import MPI
from gprMax.utilities.mpi import Dir
def calculate_starts_and_subshape(
shape: npt.NDArray[np.int32],
negative_offset: npt.NDArray[np.int32],
positive_offset: npt.NDArray[np.int32],
dirs: npt.NDArray[np.int32],
sending: bool = False,
) -> Tuple[npt.NDArray[np.int32], npt.NDArray[np.int32]]:
negative_offset = np.where(
dirs == Dir.NONE,
np.maximum(negative_offset, 0),
np.abs(negative_offset),
)
positive_offset = np.where(
dirs == Dir.NONE,
np.maximum(positive_offset, 0),
np.abs(positive_offset),
)
starts = np.select(
[dirs == Dir.NEG, dirs == Dir.POS],
[0, shape - positive_offset - sending],
default=negative_offset,
)
subshape = np.select(
[dirs == Dir.NEG, dirs == Dir.POS],
[negative_offset + sending, positive_offset + sending],
default=shape - negative_offset - positive_offset,
)
return starts, subshape
def create_mpi_type(
shape: npt.NDArray[np.int32],
negative_offset: npt.NDArray[np.int32],
positive_offset: npt.NDArray[np.int32],
dirs: npt.NDArray[np.int32],
sending: bool = False,
) -> MPI.Datatype:
starts, subshape = calculate_starts_and_subshape(
shape, negative_offset, positive_offset, dirs, sending
)
print(
f"create_mpi_type: shape={shape}, negative_offset={negative_offset}, positive_offset={positive_offset}, dirs={dirs}, sending={sending}"
)
print(f"create_mpi_type: starts={starts}, subshape={subshape}")
mpi_type = MPI.FLOAT.Create_subarray(shape.tolist(), subshape.tolist(), starts.tolist())
mpi_type.Commit()
return mpi_type

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

@@ -21,7 +21,7 @@ import itertools
import logging
import sys
from collections import OrderedDict
from typing import Any, Iterable, List, Tuple, Union
from typing import Any, Iterable, List, Optional, Tuple, Union
import numpy as np
import numpy.typing as npt
@@ -32,7 +32,8 @@ from typing_extensions import TypeVar
from gprMax import config
from gprMax.cython.pml_build import pml_average_er_mr
from gprMax.cython.yee_cell_build import build_electric_components, build_magnetic_components
from gprMax.fractals import FractalVolume
from gprMax.fractals.fractal_surface import FractalSurface
from gprMax.fractals.fractal_volume import FractalVolume
from gprMax.materials import ListMaterial, Material, PeplinskiSoil, RangeMaterial, process_materials
from gprMax.pml import CFS, PML, print_pml_info
from gprMax.receivers import Rx
@@ -185,6 +186,34 @@ class FDTDGrid:
self.pmls["thickness"]["ymax"] = int(thickness[4])
self.pmls["thickness"]["zmax"] = int(thickness[5])
def add_fractal_volume(
self,
xs: int,
xf: int,
ys: int,
yf: int,
zs: int,
zf: int,
frac_dim: float,
seed: Optional[int],
) -> FractalVolume:
volume = FractalVolume(xs, xf, ys, yf, zs, zf, frac_dim, seed)
self.fractalvolumes.append(volume)
return volume
def create_fractal_surface(
self,
xs: int,
xf: int,
ys: int,
yf: int,
zs: int,
zf: int,
frac_dim: float,
seed: Optional[int],
) -> FractalSurface:
return FractalSurface(xs, xf, ys, yf, zs, zf, frac_dim, seed)
def add_source(self, source: Source):
if isinstance(source, VoltageSource):
self.voltagesources.append(source)

30
gprMax/grid/mpi_grid.py
查看文件

@@ -28,6 +28,8 @@ from numpy import ndarray
from gprMax import config
from gprMax.cython.pml_build import pml_sum_er_mr
from gprMax.fractals.fractal_surface import MPIFractalSurface
from gprMax.fractals.fractal_volume import MPIFractalVolume
from gprMax.grid.fdtd_grid import FDTDGrid
from gprMax.pml import MPIPML, PML
from gprMax.receivers import Rx
@@ -130,6 +132,34 @@ class MPIGrid(FDTDGrid):
if self.has_neighbour(Dim.Z, Dir.POS):
self.pmls["thickness"]["zmax"] = 0
def add_fractal_volume(
self,
xs: int,
xf: int,
ys: int,
yf: int,
zs: int,
zf: int,
frac_dim: float,
seed: Optional[int],
) -> MPIFractalVolume:
volume = MPIFractalVolume(xs, xf, ys, yf, zs, zf, frac_dim, seed, self.comm, self.size)
self.fractalvolumes.append(volume)
return volume
def create_fractal_surface(
self,
xs: int,
xf: int,
ys: int,
yf: int,
zs: int,
zf: int,
frac_dim: float,
seed: Optional[int],
) -> MPIFractalSurface:
return MPIFractalSurface(xs, xf, ys, yf, zs, zf, frac_dim, seed, self.comm, self.size)
def is_coordinator(self) -> bool:
"""Test if the current rank is the coordinator.

12
gprMax/hash_cmds_geometry.py
查看文件

@@ -421,6 +421,10 @@ def process_geometrycmds(geometry):
)
raise ValueError
# Only build objects attached to the current fractal box
if tmp[12] != ID:
continue
p1 = (float(tmp[1]), float(tmp[2]), float(tmp[3]))
p2 = (float(tmp[4]), float(tmp[5]), float(tmp[6]))
frac_dim = float(tmp[7])
@@ -462,6 +466,10 @@ def process_geometrycmds(geometry):
)
raise ValueError
# Only build objects attached to the current fractal box
if tmp[8] != ID:
continue
p1 = (float(tmp[1]), float(tmp[2]), float(tmp[3]))
p2 = (float(tmp[4]), float(tmp[5]), float(tmp[6]))
depth = float(tmp[7])
@@ -477,6 +485,10 @@ def process_geometrycmds(geometry):
)
raise ValueError
# Only build objects attached to the current fractal box
if tmp[11] != ID:
continue
p1 = (float(tmp[1]), float(tmp[2]), float(tmp[3]))
p2 = (float(tmp[4]), float(tmp[5]), float(tmp[6]))
frac_dim = float(tmp[7])

129
gprMax/user_objects/cmds_geometry/add_grass.py
查看文件

@@ -20,7 +20,8 @@ import logging
import numpy as np
from gprMax.fractals import FractalSurface, Grass
from gprMax.fractals.fractal_surface import FractalSurface
from gprMax.fractals.grass import Grass
from gprMax.grid.fdtd_grid import FDTDGrid
from gprMax.materials import create_grass
from gprMax.user_objects.rotatable import RotatableMixin
@@ -92,119 +93,84 @@ class AddGrass(RotatableMixin, GeometryUserObject):
# Get the correct fractal volume
volumes = [volume for volume in grid.fractalvolumes if volume.ID == fractal_box_id]
try:
if volumes:
volume = volumes[0]
except NameError:
logger.exception(f"{self.__str__()} cannot find FractalBox {fractal_box_id}")
raise
else:
raise ValueError(f"{self.__str__()} cannot find FractalBox {fractal_box_id}")
uip = self._create_uip(grid)
_, p1, p2 = uip.check_box_points(p1, p2, self.__str__())
xs, ys, zs = p1
xf, yf, zf = p2
discretised_p1, discretised_p2 = uip.check_output_object_bounds(p1, p2, self.__str__())
xs, ys, zs = discretised_p1
xf, yf, zf = discretised_p2
if frac_dim < 0:
logger.exception(
raise ValueError(
f"{self.__str__()} requires a positive value for the fractal dimension"
)
raise ValueError
if limits[0] < 0 or limits[1] < 0:
logger.exception(
raise ValueError(
f"{self.__str__()} requires a positive value for the minimum and maximum heights for grass blades"
)
raise ValueError
# Check for valid orientations
if np.count_nonzero(discretised_p1 == discretised_p2) != 1:
raise ValueError(f"{self.__str__()} dimensions are not specified correctly")
if xs == xf:
if ys == yf or zs == zf:
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"
)
raise ValueError
fractalrange = (
round_value(limits[0] / grid.dx),
round_value(limits[1] / grid.dx),
)
# xminus surface
if xs == volume.xs:
logger.exception(
raise ValueError(
f"{self.__str__()} grass can only be specified on surfaces in the positive axis direction"
)
raise ValueError
# xplus surface
elif xf == volume.xf:
if fractalrange[1] > grid.nx:
logger.exception(
f"{self.__str__()} cannot apply grass to "
"fractal box as it would exceed the domain "
"size in the x direction"
)
raise ValueError
lower_bound = uip.discretise_point((limits[0], 0, 0))
upper_bound = uip.discretise_point((limits[1], p2[1], p2[2]))
uip.point_within_bounds(upper_bound, self.__str__())
fractalrange = (lower_bound[0], upper_bound[0])
requestedsurface = "xplus"
else:
raise ValueError(
f"{self.__str__()} must specify external surfaces on a fractal box"
)
elif ys == yf:
if zs == zf:
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"
)
raise ValueError
fractalrange = (
round_value(limits[0] / grid.dy),
round_value(limits[1] / grid.dy),
)
# yminus surface
if ys == volume.ys:
logger.exception(
raise ValueError(
f"{self.__str__()} grass can only be specified on surfaces in the positive axis direction"
)
raise ValueError
# yplus surface
elif yf == volume.yf:
if fractalrange[1] > grid.ny:
logger.exception(
f"{self.__str__()} cannot apply grass to "
"fractal box as it would exceed the domain "
"size in the y direction"
)
raise ValueError
lower_bound = uip.discretise_point((0, limits[0], 0))
upper_bound = uip.discretise_point((p2[0], limits[1], p2[2]))
uip.point_within_bounds(upper_bound, self.__str__())
fractalrange = (lower_bound[1], upper_bound[1])
requestedsurface = "yplus"
elif zs == zf:
if zs not in [volume.zs, volume.zf]:
logger.exception(
else:
raise ValueError(
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),
)
elif zs == zf:
# zminus surface
if zs == volume.zs:
logger.exception(
raise ValueError(
f"{self.__str__()} grass can only be specified on surfaces in the positive axis direction"
)
raise ValueError
# zplus surface
elif zf == volume.zf:
if fractalrange[1] > grid.nz:
logger.exception(
f"{self.__str__()} cannot apply grass to "
"fractal box as it would exceed the domain "
"size in the z direction"
)
raise ValueError
lower_bound = uip.discretise_point((0, 0, limits[0]))
upper_bound = uip.discretise_point((p2[0], p2[1], limits[1]))
uip.point_within_bounds(upper_bound, self.__str__())
fractalrange = (lower_bound[2], upper_bound[2])
requestedsurface = "zplus"
else:
raise ValueError(
f"{self.__str__()} must specify external surfaces on a fractal box"
)
else:
logger.exception(f"{self.__str__()} dimensions are not specified correctly")
raise ValueError
raise ValueError(f"{self.__str__()} dimensions are not specified correctly")
surface = FractalSurface(xs, xf, ys, yf, zs, zf, frac_dim, seed)
surface.ID = "grass"
@@ -215,11 +181,10 @@ class AddGrass(RotatableMixin, GeometryUserObject):
surface.operatingonID = volume.ID
surface.generate_fractal_surface()
if n_blades > surface.fractalsurface.shape[0] * surface.fractalsurface.shape[1]:
logger.exception(
raise ValueError(
f"{self.__str__()} the specified surface is not large "
"enough for the number of grass blades/roots specified"
)
raise ValueError
# Scale the distribution so that the summation is equal to one,
# i.e. a probability distribution
@@ -266,18 +231,20 @@ class AddGrass(RotatableMixin, GeometryUserObject):
grass = next((x for x in grid.materials if x.ID == "grass"))
testgrass = next((x for x in grass.tau if x < grid.dt), None)
if testgrass:
logger.exception(
raise ValueError(
f"{self.__str__()} requires the time step for the "
"model to be less than the relaxation time required to model grass."
)
raise ValueError
volume.fractalsurfaces.append(surface)
p3 = uip.round_to_grid_static_point(p1)
p4 = uip.round_to_grid_static_point(p2)
logger.info(
f"{self.grid_name(grid)}{n_blades} blades of grass on surface from "
f"{xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, "
f"to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m "
f"{p3[0]:g}m, {p3[1]:g}m, {p3[2]:g}m, "
f"to {p4[0]:g}m, {p4[1]:g}m, {p4[2]:g}m "
f"with fractal dimension {surface.dimension:g}, fractal seeding "
f"{surface.seed}, and range {limits[0]:g}m to {limits[1]:g}m, "
f"added to {surface.operatingonID}."

214
gprMax/user_objects/cmds_geometry/add_surface_roughness.py
查看文件

@@ -20,11 +20,9 @@ import logging
import numpy as np
from gprMax.fractals import FractalSurface
from gprMax.grid.fdtd_grid import FDTDGrid
from gprMax.user_objects.rotatable import RotatableMixin
from gprMax.user_objects.user_objects import GeometryUserObject
from gprMax.utilities.utilities import round_value
from .cmds_geometry import rotate_2point_object
@@ -37,18 +35,19 @@ class AddSurfaceRoughness(RotatableMixin, GeometryUserObject):
Attributes:
p1: list of the lower left (x,y,z) coordinates of a surface on a
FractalBox class.
p2: list of the upper right (x,y,z) coordinates of a surface on a
FractalBox class.
frac_dim: float for the fractal dimension which, for an orthogonal
parallelepiped, should take values between zero and three.
weighting: list with weightings in the first and second direction of
the surface.
limits: ist to define lower and upper limits for a range over which
the surface roughness can vary.
fractal_box_id: string identifier for the FractalBox class
that the surface roughness should be applied to.
seed: (optional) float parameter which controls the seeding of the random
number generator used to create the fractals.
p2: list of the upper right (x,y,z) coordinates of a surface on
a FractalBox class.
frac_dim: float for the fractal dimension which, for an
orthogonal parallelepiped, should take values between zero
and three.
weighting: list with weightings in the first and second
direction of the surface.
limits: list to define lower and upper limits for a range over
which the surface roughness can vary.
fractal_box_id: string identifier for the FractalBox class that
the surface roughness should be applied to.
seed: (optional) float parameter which controls the seeding of
the random number generator used to create the fractals.
"""
@property
@@ -95,145 +94,111 @@ class AddSurfaceRoughness(RotatableMixin, GeometryUserObject):
if volumes:
volume = volumes[0]
else:
logger.exception(f"{self.__str__()} cannot find FractalBox {fractal_box_id}")
raise ValueError
raise ValueError(f"{self.__str__()} cannot find FractalBox {fractal_box_id}")
uip = self._create_uip(grid)
_, p1, p2 = uip.check_box_points(p1, p2, self.__str__())
xs, ys, zs = p1
xf, yf, zf = p2
discretised_p1, discretised_p2 = uip.check_output_object_bounds(p1, p2, self.__str__())
xs, ys, zs = discretised_p1
xf, yf, zf = discretised_p2
if frac_dim < 0:
logger.exception(
raise ValueError(
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 first direction of the surface"
raise ValueError(
f"{self.__str__()} requires a positive value for the fractal weighting in the first"
" direction of the surface"
)
raise ValueError
if weighting[1] < 0:
logger.exception(
f"{self.__str__()} requires a positive value for the "
"fractal weighting in the second direction of the surface"
raise ValueError(
f"{self.__str__()} requires a positive value for the fractal weighting in the"
" second direction of the surface"
)
raise ValueError
# Check for valid orientations
if np.count_nonzero(discretised_p1 == discretised_p2) != 1:
raise ValueError(f"{self.__str__()} dimensions are not specified correctly")
if xs == xf:
if ys == yf or zs == zf:
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"
)
raise ValueError
fractalrange = (
round_value(limits[0] / grid.dx),
round_value(limits[1] / grid.dx),
)
# xminus surface
if xs == volume.xs:
if fractalrange[0] < 0 or fractalrange[1] > volume.xf:
logger.exception(
f"{self.__str__()} cannot apply fractal surface "
"to fractal box as it would exceed either the "
"upper coordinates of the fractal box or the "
"domain in the x direction"
)
raise ValueError
lower_bound = discretised_p1
upper_bound = uip.discretise_point((limits[1], p2[1], p2[2]))
grid_bound = uip.discretise_point((limits[0], p1[1], p1[2]))
fractalrange = (grid_bound[0], upper_bound[0])
requestedsurface = "xminus"
# xplus surface
elif xf == volume.xf:
if fractalrange[0] < volume.xs or fractalrange[1] > grid.nx:
logger.exception(
f"{self.__str__()} cannot apply fractal surface "
"to fractal box as it would exceed either the "
"lower coordinates of the fractal box or the "
"domain in the x direction"
)
raise ValueError
lower_bound = uip.discretise_point((limits[0], p1[1], p1[2]))
upper_bound = discretised_p2
grid_bound = uip.discretise_point((limits[1], p2[1], p2[2]))
fractalrange = (lower_bound[0], grid_bound[0])
requestedsurface = "xplus"
elif ys == yf:
if zs == zf:
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"
else:
raise ValueError(
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),
)
elif ys == yf:
# yminus surface
if ys == volume.ys:
if fractalrange[0] < 0 or fractalrange[1] > volume.yf:
logger.exception(
f"{self.__str__()} cannot apply fractal surface "
"to fractal box as it would exceed either the "
"upper coordinates of the fractal box or the "
"domain in the y direction"
)
raise ValueError
lower_bound = discretised_p1
upper_bound = uip.discretise_point((p2[0], limits[1], p2[2]))
grid_bound = uip.discretise_point((p1[0], limits[0], p1[2]))
fractalrange = (grid_bound[1], upper_bound[1])
requestedsurface = "yminus"
# yplus surface
elif yf == volume.yf:
if fractalrange[0] < volume.ys or fractalrange[1] > grid.ny:
logger.exception(
f"{self.__str__()} cannot apply fractal surface "
"to fractal box as it would exceed either the "
"lower coordinates of the fractal box or the "
"domain in the y direction"
)
raise ValueError
lower_bound = uip.discretise_point((p1[0], limits[0], p1[2]))
upper_bound = discretised_p2
grid_bound = uip.discretise_point((p2[0], limits[1], p2[2]))
fractalrange = (lower_bound[1], grid_bound[1])
requestedsurface = "yplus"
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"
else:
raise ValueError(
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),
)
elif zs == zf:
# zminus surface
if zs == volume.zs:
if fractalrange[0] < 0 or fractalrange[1] > volume.zf:
logger.exception(
f"{self.__str__()} cannot apply fractal surface "
"to fractal box as it would exceed either the "
"upper coordinates of the fractal box or the "
"domain in the x direction"
)
raise ValueError
lower_bound = discretised_p1
upper_bound = uip.discretise_point((p2[0], p2[1], limits[1]))
grid_bound = uip.discretise_point((p1[0], p1[1], limits[0]))
fractalrange = (grid_bound[2], upper_bound[2])
requestedsurface = "zminus"
# zplus surface
elif zf == volume.zf:
if fractalrange[0] < volume.zs or fractalrange[1] > grid.nz:
logger.exception(
f"{self.__str__()} cannot apply fractal surface "
"to fractal box as it would exceed either the "
"lower coordinates of the fractal box or the "
"domain in the z direction"
)
raise ValueError
lower_bound = uip.discretise_point((p1[0], p1[1], limits[0]))
upper_bound = discretised_p2
grid_bound = uip.discretise_point((p2[0], p2[1], limits[1]))
fractalrange = (lower_bound[2], grid_bound[2])
requestedsurface = "zplus"
else:
raise ValueError(
f"{self.__str__()} can only be used on the external surfaces of a fractal box"
)
else:
logger.exception(f"{self.__str__()} dimensions are not specified correctly")
raise ValueError
raise ValueError(f"{self.__str__()} dimensions are not specified correctly")
surface = FractalSurface(xs, xf, ys, yf, zs, zf, frac_dim, seed)
if any(lower_bound < volume.start):
raise ValueError(
f"{self.__str__()} cannot apply fractal surface to"
" fractal box as it would exceed the lower coordinates"
" of the fractal box."
)
elif any(upper_bound > volume.stop):
raise ValueError(
f"{self.__str__()} cannot apply fractal surface to"
" fractal box as it would exceed the upper coordinates"
" of the fractal box."
)
# Check lower or upper extent of the fractal surface (depending
# if the fractal surface has been applied in the minus or plus
# direction).
uip.point_within_bounds(grid_bound, f"{self.__str__()}")
surface = grid.create_fractal_surface(xs, xf, ys, yf, zs, zf, frac_dim, seed)
surface.surfaceID = requestedsurface
surface.fractalrange = fractalrange
surface.operatingonID = volume.ID
@@ -242,18 +207,21 @@ class AddSurfaceRoughness(RotatableMixin, GeometryUserObject):
# List of existing surfaces IDs
existingsurfaceIDs = [x.surfaceID for x in volume.fractalsurfaces]
if surface.surfaceID in existingsurfaceIDs:
logger.exception(
raise ValueError(
f"{self.__str__()} has already been used on the {surface.surfaceID} surface"
)
raise ValueError
surface.generate_fractal_surface()
volume.fractalsurfaces.append(surface)
p3 = uip.round_to_grid_static_point(p1)
p4 = uip.round_to_grid_static_point(p2)
logger.info(
f"{self.grid_name(grid)}Fractal surface from {xs * grid.dx:g}m, "
f"{ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, "
f"{yf * grid.dy:g}m, {zf * grid.dz:g}m with fractal dimension "
f"{self.grid_name(grid)}Fractal surface from {p3[0]:g}m, "
f"{p3[1]:g}m, {p3[2]:g}m, to {p4[0]:g}m, "
f"{p4[1]:g}m, {p4[2]:g}m with fractal dimension "
f"{surface.dimension:g}, fractal weightings {surface.weighting[0]:g}, "
f"{surface.weighting[1]:g}, fractal seeding {surface.seed}, "
f"and range {limits[0]:g}m to {limits[1]:g}m, added to "

84
gprMax/user_objects/cmds_geometry/add_surface_water.py
查看文件

@@ -24,7 +24,6 @@ from gprMax.grid.fdtd_grid import FDTDGrid
from gprMax.materials import create_water
from gprMax.user_objects.rotatable import RotatableMixin
from gprMax.user_objects.user_objects import GeometryUserObject
from gprMax.utilities.utilities import round_value
from .cmds_geometry import rotate_2point_object
@@ -77,80 +76,70 @@ class AddSurfaceWater(RotatableMixin, GeometryUserObject):
if volumes := [volume for volume in grid.fractalvolumes if volume.ID == fractal_box_id]:
volume = volumes[0]
else:
logger.exception(f"{self.__str__()} cannot find FractalBox {fractal_box_id}")
raise ValueError
raise ValueError(f"{self.__str__()} cannot find FractalBox {fractal_box_id}")
uip = self._create_uip(grid)
_, p1, p2 = uip.check_box_points(p1, p2, self.__str__())
xs, ys, zs = p1
xf, yf, zf = p2
discretised_p1, discretised_p2 = uip.check_output_object_bounds(p1, p2, self.__str__())
xs, ys, zs = discretised_p1
xf, yf, zf = discretised_p2
if depth <= 0:
logger.exception(f"{self.__str__()} requires a positive value for the depth of water")
raise ValueError
raise ValueError(f"{self.__str__()} requires a positive value for the depth of water")
# Check for valid orientations
if np.count_nonzero(discretised_p1 == discretised_p2) != 1:
raise ValueError(f"{self.__str__()} dimensions are not specified correctly")
if xs == xf:
if ys == yf or zs == zf:
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"
)
raise ValueError
# xminus surface
if xs == volume.xs:
requestedsurface = "xminus"
# xplus surface
elif xf == volume.xf:
requestedsurface = "xplus"
filldepthcells = round_value(depth / grid.dx)
filldepth = filldepthcells * grid.dx
elif ys == yf:
if zs == zf:
logger.exception(f"{self.__str__()} dimensions are not specified correctly")
raise ValueError
if ys not in [volume.ys, volume.yf]:
logger.exception(
else:
raise ValueError(
f"{self.__str__()} can only be used on the external surfaces of a fractal box"
)
raise ValueError
filldepthcells = uip.discretise_point((depth, 0, 0))[0]
filldepth = uip.round_to_grid_static_point((depth, 0, 0))[0]
elif ys == yf:
# yminus surface
if ys == volume.ys:
requestedsurface = "yminus"
# yplus surface
elif yf == volume.yf:
requestedsurface = "yplus"
filldepthcells = round_value(depth / grid.dy)
filldepth = filldepthcells * grid.dy
elif zs == zf:
if zs not in [volume.zs, volume.zf]:
logger.exception(
else:
raise ValueError(
f"{self.__str__()} can only be used on the external surfaces of a fractal box"
)
raise ValueError
filldepthcells = uip.discretise_point((0, depth, 0))[1]
filldepth = uip.round_to_grid_static_point((0, depth, 0))[1]
elif zs == zf:
# zminus surface
if zs == volume.zs:
requestedsurface = "zminus"
# zplus surface
elif zf == volume.zf:
requestedsurface = "zplus"
filldepthcells = round_value(depth / grid.dz)
filldepth = filldepthcells * grid.dz
else:
raise ValueError(
f"{self.__str__()} can only be used on the external surfaces of a fractal box"
)
filldepthcells = uip.discretise_point((0, 0, depth))[2]
filldepth = uip.round_to_grid_static_point((0, 0, depth))[2]
else:
logger.exception(f"{self.__str__()} dimensions are not specified correctly")
raise ValueError
raise ValueError(f"{self.__str__()} dimensions are not specified correctly")
surface = next((x for x in volume.fractalsurfaces if x.surfaceID == requestedsurface), None)
if not surface:
logger.exception(
raise ValueError(
f"{self.__str__()} specified surface {requestedsurface} does not have a rough surface applied"
)
raise ValueError
surface.filldepth = filldepthcells
@@ -159,11 +148,10 @@ class AddSurfaceWater(RotatableMixin, GeometryUserObject):
surface.filldepth < surface.fractalrange[0]
or surface.filldepth > surface.fractalrange[1]
):
logger.exception(
raise ValueError(
f"{self.__str__()} requires a value for the depth of water that lies with the "
f"range of the requested surface roughness"
)
raise ValueError
# Check to see if water has been already defined as a material
if not any(x.ID == "water" for x in grid.materials):
@@ -172,15 +160,17 @@ class AddSurfaceWater(RotatableMixin, GeometryUserObject):
# Check if time step for model is suitable for using water
water = next((x for x in grid.materials if x.ID == "water"))
if testwater := next((x for x in water.tau if x < grid.dt), None):
logger.exception(
raise ValueError(
f"{self.__str__()} requires the time step for the model "
f"to be less than the relaxation time required to model water."
)
raise ValueError
p3 = uip.round_to_grid_static_point(p1)
p4 = uip.round_to_grid_static_point(p2)
logger.info(
f"{self.grid_name(grid)}Water on surface from {xs * grid.dx:g}m, "
f"{ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, "
f"{yf * grid.dy:g}m, {zf * grid.dz:g}m with depth {filldepth:g}m, "
f"added to {surface.operatingonID}."
f"{self.grid_name(grid)}Water on surface from {p3[0]:g}m,"
f" {p3[1]:g}m, {p3[2]:g}m, to {p4[0]:g}m, {p4[1]:g}m,"
f" {p4[2]:g}m with depth {filldepth:g}m, added to"
f" {surface.operatingonID}."
)

53
gprMax/user_objects/cmds_geometry/fractal_box.py
查看文件

@@ -22,7 +22,6 @@ import numpy as np
import gprMax.config as config
from gprMax.cython.geometry_primitives import build_voxels_from_array, build_voxels_from_array_mask
from gprMax.fractals import FractalVolume
from gprMax.grid.fdtd_grid import FDTDGrid
from gprMax.materials import ListMaterial
from gprMax.user_objects.cmds_geometry.cmds_geometry import check_averaging, rotate_2point_object
@@ -110,12 +109,7 @@ class FractalBox(RotatableMixin, GeometryUserObject):
p3 = uip.round_to_grid_static_point(p1)
p4 = uip.round_to_grid_static_point(p2)
grid_contains_fractal_box, p1, p2 = uip.check_box_points(p1, p2, self.__str__())
# Exit early if none of the fractal box is in this grid as there
# is nothing else to do.
if not grid_contains_fractal_box:
return
p1, p2 = uip.check_output_object_bounds(p1, p2, self.__str__())
xs, ys, zs = p1
xf, yf, zf = p2
@@ -172,7 +166,7 @@ class FractalBox(RotatableMixin, GeometryUserObject):
)
raise ValueError
self.volume = FractalVolume(xs, xf, ys, yf, zs, zf, frac_dim, seed)
self.volume = grid.add_fractal_volume(xs, xf, ys, yf, zs, zf, frac_dim, seed)
self.volume.ID = ID
self.volume.operatingonID = mixing_model_id
self.volume.nbins = nbins
@@ -199,39 +193,26 @@ class FractalBox(RotatableMixin, GeometryUserObject):
self.do_pre_build = False
else:
if self.volume.fractalsurfaces:
self.volume.originalxs = self.volume.xs
self.volume.originalxf = self.volume.xf
self.volume.originalys = self.volume.ys
self.volume.originalyf = self.volume.yf
self.volume.originalzs = self.volume.zs
self.volume.originalzf = self.volume.zf
# Extend the volume to accomodate any rough surfaces, grass,
# or roots
for surface in self.volume.fractalsurfaces:
if surface.surfaceID == "xminus":
if surface.fractalrange[0] < self.volume.xs:
self.volume.nx += self.volume.xs - surface.fractalrange[0]
self.volume.xs = surface.fractalrange[0]
elif surface.surfaceID == "xplus":
if surface.fractalrange[1] > self.volume.xf:
self.volume.nx += surface.fractalrange[1] - self.volume.xf
self.volume.xf = surface.fractalrange[1]
elif surface.surfaceID == "yminus":
if surface.fractalrange[0] < self.volume.ys:
self.volume.ny += self.volume.ys - surface.fractalrange[0]
self.volume.ys = surface.fractalrange[0]
elif surface.surfaceID == "yplus":
if surface.fractalrange[1] > self.volume.yf:
self.volume.ny += surface.fractalrange[1] - self.volume.yf
self.volume.yf = surface.fractalrange[1]
elif surface.surfaceID == "zminus":
if surface.fractalrange[0] < self.volume.zs:
self.volume.nz += self.volume.zs - surface.fractalrange[0]
self.volume.zs = surface.fractalrange[0]
elif surface.surfaceID == "zplus":
if surface.fractalrange[1] > self.volume.zf:
self.volume.nz += surface.fractalrange[1] - self.volume.zf
self.volume.zf = surface.fractalrange[1]
# If there is only 1 bin then a normal material is being used,
@@ -246,7 +227,8 @@ class FractalBox(RotatableMixin, GeometryUserObject):
)
self.volume.fractalvolume *= materialnumID
else:
self.volume.generate_fractal_volume()
if not self.volume.generate_fractal_volume():
return
for i in range(0, self.volume.nx):
for j in range(0, self.volume.ny):
for k in range(0, self.volume.nz):
@@ -262,6 +244,10 @@ class FractalBox(RotatableMixin, GeometryUserObject):
# TODO: Allow extract of rough surface profile (to print/file?)
for surface in self.volume.fractalsurfaces:
if surface.surfaceID == "xminus":
surface.fractalrange = (
max(surface.fractalrange[0], 0),
min(surface.fractalrange[1], grid.nx),
)
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):
@@ -285,6 +271,10 @@ class FractalBox(RotatableMixin, GeometryUserObject):
] = 0
elif surface.surfaceID == "xplus":
surface.fractalrange = (
max(surface.fractalrange[0], 0),
min(surface.fractalrange[1], grid.nx),
)
if not surface.ID:
for i in range(surface.fractalrange[0], surface.fractalrange[1]):
for j in range(surface.ys, surface.yf):
@@ -393,6 +383,10 @@ class FractalBox(RotatableMixin, GeometryUserObject):
root += 1
elif surface.surfaceID == "yminus":
surface.fractalrange = (
max(surface.fractalrange[0], 0),
min(surface.fractalrange[1], grid.ny),
)
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):
@@ -416,6 +410,10 @@ class FractalBox(RotatableMixin, GeometryUserObject):
] = 0
elif surface.surfaceID == "yplus":
surface.fractalrange = (
max(surface.fractalrange[0], 0),
min(surface.fractalrange[1], grid.ny),
)
if not surface.ID:
for i in range(surface.xs, surface.xf):
for j in range(surface.fractalrange[0], surface.fractalrange[1]):
@@ -524,6 +522,10 @@ class FractalBox(RotatableMixin, GeometryUserObject):
root += 1
elif surface.surfaceID == "zminus":
surface.fractalrange = (
max(surface.fractalrange[0], 0),
min(surface.fractalrange[1], grid.nz),
)
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]):
@@ -547,6 +549,10 @@ class FractalBox(RotatableMixin, GeometryUserObject):
] = 0
elif surface.surfaceID == "zplus":
surface.fractalrange = (
max(surface.fractalrange[0], 0),
min(surface.fractalrange[1], grid.nz),
)
if not surface.ID:
for i in range(surface.xs, surface.xf):
for j in range(surface.ys, surface.yf):
@@ -685,7 +691,8 @@ class FractalBox(RotatableMixin, GeometryUserObject):
)
raise ValueError
else:
self.volume.generate_fractal_volume()
if not self.volume.generate_fractal_volume():
return
for i in range(0, self.volume.nx):
for j in range(0, self.volume.ny):
for k in range(0, self.volume.nz):

47
gprMax/utilities/mpi.py 普通文件
查看文件

@@ -0,0 +1,47 @@
from enum import IntEnum, unique
from typing import Union
import numpy as np
import numpy.typing as npt
from mpi4py import MPI
@unique
class Dim(IntEnum):
X = 0
Y = 1
Z = 2
@unique
class Dir(IntEnum):
NONE = -1
NEG = 0
POS = 1
def get_neighbours(comm: MPI.Cartcomm) -> npt.NDArray[np.int32]:
neighbours = np.full((3, 2), -1, dtype=np.int32)
neighbours[Dim.X] = comm.Shift(direction=Dim.X, disp=1)
neighbours[Dim.Y] = comm.Shift(direction=Dim.Y, disp=1)
neighbours[Dim.Z] = comm.Shift(direction=Dim.Z, disp=1)
return neighbours
def get_neighbour(comm: MPI.Cartcomm, dim: Dim, dir: Dir) -> int:
neighbours = comm.Shift(direction=dim, disp=1)
return neighbours[dir]
def get_relative_neighbour(
comm: MPI.Cartcomm, dirs: npt.NDArray[np.int32], disp: Union[int, npt.NDArray[np.int32]] = 1
) -> int:
offset = np.select([dirs == Dir.NEG, dirs == Dir.POS], [-disp, disp], default=0)
coord = comm.coords + offset
if any(coord < 0) or any(coord >= comm.dims):
return -1
return comm.Get_cart_rank(coord.tolist())

二进制
reframe_tests/regression_checks/TestAddSurfaceRoughness_2a49ee92/add_surface_roughness_boundaries.h5 普通文件
查看文件

二进制文件未显示。

4
reframe_tests/regression_checks/TestFractalBoxGeometry_1887bd81/full_volume_materials.txt → reframe_tests/regression_checks/TestAddSurfaceRoughness_2a49ee92/full_volume_materials.txt
查看文件

@@ -11,12 +11,12 @@
#add_dispersion_debye: 1 8.6327 8.0994e-12 |10.5897+0.0769+1.0000+0.0000|
#material: 10.7105 0.0777078 1 0 |10.7105+0.0777+1.0000+0.0000|
#add_dispersion_debye: 1 8.9052 8.0994e-12 |10.7105+0.0777+1.0000+0.0000|
#material: 10.8308 0.0784642 1 0 |10.8308+0.0785+1.0000+0.0000|
#add_dispersion_debye: 1 9.18039 8.0994e-12 |10.8308+0.0785+1.0000+0.0000|
#material: 10.9505 0.0792123 1 0 |10.9505+0.0792+1.0000+0.0000|
#add_dispersion_debye: 1 9.45822 8.0994e-12 |10.9505+0.0792+1.0000+0.0000|
#material: 3.67252 0.0108321 1 0 |3.6725+0.0108+1.0000+0.0000|
#add_dispersion_debye: 1 0.0182374 8.0994e-12 |3.6725+0.0108+1.0000+0.0000|
#material: 3.91944 0.0163865 1 0 |3.9194+0.0164+1.0000+0.0000|
#add_dispersion_debye: 1 0.0669568 8.0994e-12 |3.9194+0.0164+1.0000+0.0000|
#material: 4.1448 0.0203355 1 0 |4.1448+0.0203+1.0000+0.0000|
#add_dispersion_debye: 1 0.131947 8.0994e-12 |4.1448+0.0203+1.0000+0.0000|
#material: 4.35657 0.0235572 1 0 |4.3566+0.0236+1.0000+0.0000|

二进制
reframe_tests/regression_checks/TestAddSurfaceRoughness_9b260413/add_surface_roughness_half.h5
查看文件

二进制文件未显示。

二进制
reframe_tests/regression_checks/TestFractalBoxGeometry_090b493e/fractal_box_boundaries.h5 普通文件
查看文件

二进制文件未显示。

6
reframe_tests/regression_checks/TestAddSurfaceRoughness_9b260413/full_volume_materials.txt → reframe_tests/regression_checks/TestFractalBoxGeometry_090b493e/full_volume_materials.txt
查看文件

@@ -11,14 +11,12 @@
#add_dispersion_debye: 1 8.6327 8.0994e-12 |10.5897+0.0769+1.0000+0.0000|
#material: 10.7105 0.0777078 1 0 |10.7105+0.0777+1.0000+0.0000|
#add_dispersion_debye: 1 8.9052 8.0994e-12 |10.7105+0.0777+1.0000+0.0000|
#material: 10.8308 0.0784642 1 0 |10.8308+0.0785+1.0000+0.0000|
#add_dispersion_debye: 1 9.18039 8.0994e-12 |10.8308+0.0785+1.0000+0.0000|
#material: 10.9505 0.0792123 1 0 |10.9505+0.0792+1.0000+0.0000|
#add_dispersion_debye: 1 9.45822 8.0994e-12 |10.9505+0.0792+1.0000+0.0000|
#material: 3.67252 0.0108321 1 0 |3.6725+0.0108+1.0000+0.0000|
#add_dispersion_debye: 1 0.0182374 8.0994e-12 |3.6725+0.0108+1.0000+0.0000|
#material: 3.91944 0.0163865 1 0 |3.9194+0.0164+1.0000+0.0000|
#add_dispersion_debye: 1 0.0669568 8.0994e-12 |3.9194+0.0164+1.0000+0.0000|
#material: 4.1448 0.0203355 1 0 |4.1448+0.0203+1.0000+0.0000|
#add_dispersion_debye: 1 0.131947 8.0994e-12 |4.1448+0.0203+1.0000+0.0000|
#material: 4.35657 0.0235572 1 0 |4.3566+0.0236+1.0000+0.0000|
#add_dispersion_debye: 1 0.209447 8.0994e-12 |4.3566+0.0236+1.0000+0.0000|
#material: 4.55835 0.0263403 1 0 |4.5583+0.0263+1.0000+0.0000|

二进制
reframe_tests/regression_checks/TestFractalBoxGeometry_1887bd81/fractal_box_half.h5
查看文件

二进制文件未显示。

2
reframe_tests/tests/mixins.py
查看文件

@@ -29,7 +29,7 @@ else:
class ReceiverMixin(GprMaxMixin):
number_of_receivers = variable(int, value=-1)
@run_after("setup")
@run_after("setup", always_last=True)
def add_receiver_regression_checks(self):
reference_file = self.build_reference_filepath(self.output_file)

8
reframe_tests/tests/src/geometry_tests/fractal_box_geometry/fractal_box_half.in → reframe_tests/tests/src/geometry_tests/add_surface_roughness_geometry/add_surface_roughness_boundaries.in
查看文件

@@ -8,6 +8,12 @@
#rx: 0.080 0.080 0.080
#soil_peplinski: 0.5 0.5 2.0 2.66 0.001 0.25 my_soil
#fractal_box: 0 0 0 0.05 0.05 0.1 1.5 1 1 1 50 my_soil my_soil_box 42
#fractal_box: 0 0 0 0.049 0.05 0.051 1.5 1 1 1 50 my_soil box_1 42
#fractal_box: 0.049 0.05 0.051 0.1 0.1 0.1 1.5 1 1 1 50 my_soil box_2 22
#add_surface_roughness: 0.049 0 0 0.049 0.05 0.05 1 1 1 0.025 0.05 box_1 43
#add_surface_roughness: 0.049 0.05 0.051 0.049 0.1 0.1 1 1 1 0.049 0.074 box_2 23
#geometry_objects_write: 0 0 0 0.1 0.1 0.1 full_volume

4
reframe_tests/tests/src/geometry_tests/add_surface_roughness_geometry/add_surface_roughness_half.in → reframe_tests/tests/src/geometry_tests/fractal_box_geometry/fractal_box_boundaries.in
查看文件

@@ -8,7 +8,7 @@
#rx: 0.080 0.080 0.080
#soil_peplinski: 0.5 0.5 2.0 2.66 0.001 0.25 my_soil
#fractal_box: 0 0 0 0.05 0.05 0.1 1.5 1 1 1 50 my_soil my_soil_box 42
#add_surface_roughness: 0 0 0.1 0.05 0.05 0.1 1 1 1 0.05 0.1 my_soil_box 22
#fractal_box: 0 0 0 0.049 0.05 0.051 1.5 1 1 1 50 my_soil my_soil_box 42
#fractal_box: 0.049 0.05 0.051 0.1 0.1 0.1 1.5 1 1 1 50 my_soil my_soil_box 22
#geometry_objects_write: 0 0 0 0.1 0.1 0.1 full_volume

101
reframe_tests/tests/test_geometry.py
查看文件

@@ -8,6 +8,13 @@ from reframe_tests.tests.standard_tests import GprMaxGeometryTest, GprMaxRegress
"""
@rfm.simple_test
class TestAddGrass(GprMaxGeometryTest):
tags = {"test", "serial", "geometery", "fractal", "surface", "grass"}
sourcesdir = "src/geometry_tests/add_grass_geometry"
model = parameter(["add_grass_full", "add_grass_small"])
@rfm.simple_test
class TestAddSurfaceRoughness(GprMaxGeometryTest):
tags = {"test", "serial", "geometery", "fractal", "surface", "roughness"}
@@ -15,9 +22,9 @@ class TestAddSurfaceRoughness(GprMaxGeometryTest):
model = parameter(
[
"add_surface_roughness_full",
"add_surface_roughness_half",
"add_surface_roughness_small",
"add_surface_roughness_weighted",
"add_surface_roughness_boundaries",
]
)
@@ -29,13 +36,6 @@ class TestAddSurfaceWater(GprMaxGeometryTest):
model = parameter(["add_surface_water_full", "add_surface_water_small"])
@rfm.simple_test
class TestAddGrass(GprMaxGeometryTest):
tags = {"test", "serial", "geometery", "fractal", "surface", "grass"}
sourcesdir = "src/geometry_tests/add_grass_geometry"
model = parameter(["add_grass_full", "add_grass_small"])
@rfm.simple_test
class TestBoxGeometryDefaultPml(GprMaxRegressionTest):
tags = {"test", "serial", "geometery", "box"}
@@ -130,7 +130,12 @@ class TestFractalBoxGeometry(GprMaxGeometryTest):
tags = {"test", "serial", "geometery", "fractal", "box", "fractal_box"}
sourcesdir = "src/geometry_tests/fractal_box_geometry"
model = parameter(
["fractal_box_full", "fractal_box_half", "fractal_box_small", "fractal_box_weighted"]
[
"fractal_box_full",
"fractal_box_small",
"fractal_box_weighted",
"fractal_box_boundaries",
]
)
@@ -160,6 +165,64 @@ class TestTriangleGeometry(GprMaxRegressionTest):
"""
@rfm.simple_test
class TestAddGrassMpi(MpiMixin, TestAddGrass):
tags = {"test", "mpi", "geometery", "fractal", "surface", "grass"}
mpi_layout = parameter([[1, 2, 2], [3, 1, 3], [4, 1, 4]])
test_dependency = TestAddGrass
# Depending on the orientation/alignment of the fractal surface, not all
# MPI decompositions are valid. Therefore create tests that divide the
# domain in each pair of dimensions.
@rfm.simple_test
class TestAddSurfaceRoughnessMpiYZ(MpiMixin, TestAddSurfaceRoughness):
tags = {"test", "mpi", "geometery", "fractal", "surface", "roughness"}
mpi_layout = parameter([[1, 2, 2], [1, 3, 3], [1, 4, 4]])
test_dependency = TestAddSurfaceRoughness
model = parameter(
[
"add_surface_roughness_full",
"add_surface_roughness_small",
"add_surface_roughness_weighted",
]
)
@rfm.simple_test
class TestAddSurfaceRoughnessMpiXZ(MpiMixin, TestAddSurfaceRoughness):
tags = {"test", "mpi", "geometery", "fractal", "surface", "roughness"}
mpi_layout = parameter([[2, 1, 2], [3, 1, 3], [4, 1, 4]])
test_dependency = TestAddSurfaceRoughness
model = parameter(
[
"add_surface_roughness_full",
"add_surface_roughness_weighted",
"add_surface_roughness_boundaries",
]
)
@rfm.simple_test
class TestAddSurfaceRoughnessMpiXY(MpiMixin, TestAddSurfaceRoughness):
tags = {"test", "mpi", "geometery", "fractal", "surface", "roughness"}
mpi_layout = parameter([[2, 2, 1], [3, 3, 1], [4, 4, 1]])
test_dependency = TestAddSurfaceRoughness
model = parameter(
[
"add_surface_roughness_small",
"add_surface_roughness_boundaries",
]
)
@rfm.simple_test
class TestAddSurfaceWaterMpi(MpiMixin, TestAddSurfaceWater):
tags = {"test", "mpi", "geometery", "fractal", "surface", "roughness", "water"}
mpi_layout = parameter([[1, 2, 2], [3, 1, 3], [4, 1, 4]])
test_dependency = TestAddSurfaceWater
@rfm.simple_test
class TestBoxGeometryDefaultPmlMpi(MpiMixin, TestBoxGeometryDefaultPml):
tags = {"test", "mpi", "geometery", "box"}
@@ -219,10 +282,28 @@ class TestEllipsoidGeometryMpi(MpiMixin, TestEllipsoidGeometry):
@rfm.simple_test
class TestFractalBoxGeometryMpi(MpiMixin, TestFractalBoxGeometry):
tags = {"test", "mpi", "geometery", "fractal", "box", "fractal_box"}
mpi_layout = parameter([[2, 2, 2], [3, 3, 3], [4, 4, 4]])
mpi_layout = parameter([[2, 2, 1], [1, 3, 3], [4, 1, 4]])
test_dependency = TestFractalBoxGeometry
@rfm.simple_test
class TestFractalBoxGeometryMpiDecomposition(MpiMixin, TestFractalBoxGeometry):
"""Extra Fractal Box test.
This tests that so long as the area covered by the fractal box has
an MPI decomposition of 1 in a dimension, the model will run
successfully.
N.B: These tests would fail if run using the 'fractal_box_full'
model.
"""
tags = {"test", "mpi", "geometery", "fractal", "box", "fractal_box"}
mpi_layout = parameter([[2, 4, 4], [4, 2, 4], [4, 4, 2]])
test_dependency = TestFractalBoxGeometry
model = parameter(["fractal_box_small"])
@rfm.simple_test
class TestPlateGeometryMpi(MpiMixin, TestPlateGeometry):
tags = {"test", "mpi", "geometery", "plate"}

1
requirements.txt
查看文件

@@ -9,6 +9,7 @@ jinja2
jupyter
matplotlib
# mpi4py
mpi4py-fft
numpy<2
numpy-stl
pandas