你已经派生过 gprMax
镜像自地址
https://gitee.com/sunhf/gprMax.git
已同步 2025-08-07 15:10:13 +08:00
Merge branch 'devel-reframe' into mpi
这个提交包含在：
Nathan Mannall
2
CREDITS
2
CREDITS
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This is the official list of entities and people who have contributed to gprMax
|
||||
|
||||
8
docs/.gitignore
vendored
普通文件
8
docs/.gitignore
vendored
普通文件
@@ -0,0 +1,8 @@
|
||||
_*/
|
||||
doctrees/
|
||||
dirhtml/
|
||||
|
||||
*.html
|
||||
.buildinfo
|
||||
objects.inv
|
||||
searchindex.js
|
||||
@@ -222,6 +222,10 @@ Transmission Line
|
||||
-----------------
|
||||
.. autoclass:: gprMax.cmds_multiuse.TransmissionLine
|
||||
|
||||
Discrete Plane Wave
|
||||
-------------------
|
||||
.. autoclass:: gprMax.cmds_multiuse.DiscretePlaneWave
|
||||
|
||||
Excitation File
|
||||
---------------
|
||||
.. autoclass:: gprMax.cmds_multiuse.ExcitationFile
|
||||
|
||||
@@ -900,6 +900,29 @@ For example, to specify a z directed transmission line source with a resistance
|
||||
|
||||
An example antenna model using a transmission line can be found in the :ref:`examples <example-wire-dipole>` section.
|
||||
|
||||
#discrete_plane_wave:
|
||||
---------------------
|
||||
|
||||
Allows you to introduce a discrete plane wave source [TAN2010]_. Plane wave sources are a useful tool in multiple different scenarios of electromagnetic simulations, especially when the wave is emitted by a source that is quite far away from the target. The syntax of the command is:
|
||||
|
||||
.. code-block:: none
|
||||
|
||||
#discrete_plane_wave: f1 f2 f3 f4 f5 f6 f7 f8 f9 str1 [f10 f11]
|
||||
|
||||
* ``f1 f2 f3`` are the lower left (x,y,z) coordinates of the total field, scattered field (TFSF) box, and ``f4 f5 f6`` are the upper right (x,y,z) coordinates of the total field, scattered field (TFSF) box.
|
||||
* ``f7`` is psi which defines the polarisation of the incident plane wave.
|
||||
* ``f8`` is phi which defines the azimuthal propagation angle (degrees) of the incident plane wave.
|
||||
* ``f9`` is theta which defines the polar propagation angle (degrees) of the incident plane wave.
|
||||
* ``f10 f11`` are optional parameters. ``f10`` is a time delay in starting the excitation of the discrete plane wave. ``f11`` is a time to remove the excitation of the discrete plane wave. If the time window is longer than the excitation of the discrete plane wave removal time then the excitation of the discrete plane wave will stop after the excitation of the discrete plane wave removal time. If the excitation of the discrete plane wave removal time is longer than the time window then the excitation of the discrete plane wave will be active for the entire time window. If ``f10 f11`` are omitted the excitation of the discrete plane wave will start at the beginning of time window and stop at the end of the time window.
|
||||
* ``str1`` is the identifier of the waveform that should be used with the source.
|
||||
|
||||
For example, to specify a discrete plane wave in a TFSF box (0.010, 0.010, 0.010 to 0.040, 0.040, 0.040) with a polarisation angle of 90 degrees, azimuthal propagation angle of 63.4 degrees, polar propagation angle of 36.7 degrees, and using the waveform defined by the identifier ``mypulse`` use: ``#discrete_plane_wave: 0.010 0.010 0.010 0.040 0.040 0.040 90.0 63.4 36.7 mypulse``.
|
||||
|
||||
.. note::
|
||||
|
||||
* Currently a plane wave can be introduced from the direction of the origin only, if the Total field/scattered field box is considered to be present in the first octant (region in the coordinate space determined by the positive x-axis, positive y-axis, and positive z-axis). The implementation should ideally be able to handle the introduction of the plane wave from any direction but currently it supports 0 <= phi <= 90 degrees and 0 <= theta <= 90 degrees only.
|
||||
* The plane wave implementation was initiated through a `Google Summer of Code <https://summerofcode.withgoogle.com/>`_ (GSoC) project and `more details can be found in the original pull request <https://github.com/gprMax/gprMax/pull/373>`_.
|
||||
|
||||
#rx:
|
||||
----
|
||||
|
||||
|
||||
@@ -23,6 +23,7 @@ References
|
||||
.. [PIE2009] Pieraccini, M., Bicci, A., Mecatti, D., Macaluso, G., & Atzeni, C. (2009). Propagation of large bandwidth microwave signals in water. Antennas and Propagation, IEEE Transactions on, 57(11), 3612-3618. (http://dx.doi.org/10.1109/tap.2009.2025674)
|
||||
.. [STA2017] Stadler, S. (2017). A Forward Modeling Study for the Investigation of the Vertical Water-Content Distribution of Using Guided GPR Waves (Master’s Thesis, Freiberg University of Mining and Technology, Germany) (https://app.box.com/s/h2n2ytsdllcm77du1o9lknumngmq8vq8)
|
||||
.. [TAF2005] Taflove, A., & Hagness, S. C. (2005). Computational electrodynamics. Artech house.
|
||||
.. [TAN2010] Tan, T., & Potter, M. (2010). FDTD discrete planewave (FDTD-DPW) formulation for a perfectly matched source in TFSF simulations. IEEE Transactions on Antennas and Propagation, 58(8), 2641-2648. (https://doi.org/10.1109/TAP.2010.2050446)
|
||||
.. [TEI1998] Teixeira, F. L., Chew, W. C., Straka, M., Oristaglio, M. L., & Wang, T. (1998). Finite-difference time-domain simulation of ground penetrating radar on dispersive, inhomogeneous, and conductive soils. Geoscience and Remote Sensing, IEEE Transactions on, 36(6), 1928-1937. (http://dx.doi.org/10.1109/36.729364)
|
||||
.. [TUR1987] Turcotte, D. L. (1987). A fractal interpretation of topography and geoid spectra on the Earth, Moon, Venus, and Mars. Journal of Geophysical Research: Solid Earth (1978–2012), 92(B4), E597-E601. (http://dx.doi.org/10.1029/jb092ib04p0e597)
|
||||
.. [TUR1997] Turcotte, D. L. (1997). Fractals and chaos in geology and geophysics. Cambridge university press. (http://dx.doi.org/10.1017/cbo9781139174695)
|
||||
|
||||
@@ -32,6 +32,7 @@ from .user_objects.cmds_multiuse import (
|
||||
AddDebyeDispersion,
|
||||
AddDrudeDispersion,
|
||||
AddLorentzDispersion,
|
||||
DiscretePlaneWave,
|
||||
ExcitationFile,
|
||||
HertzianDipole,
|
||||
MagneticDipole,
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -74,7 +74,7 @@ class ModelConfig:
|
||||
deviceID = 0
|
||||
|
||||
self.device = {
|
||||
"dev": sim_config.get_model_device(deviceID),
|
||||
"dev": sim_config.set_model_device(deviceID),
|
||||
"deviceID": deviceID,
|
||||
"snapsgpu2cpu": False,
|
||||
}
|
||||
@@ -106,7 +106,11 @@ class ModelConfig:
|
||||
# phase-velocity phase error.
|
||||
# mingridsampling: minimum grid sampling of smallest wavelength for
|
||||
# physical wave propagation.
|
||||
self.numdispersion = {"highestfreqthres": 40, "maxnumericaldisp": 2, "mingridsampling": 3}
|
||||
self.numdispersion = {
|
||||
"highestfreqthres": 40,
|
||||
"maxnumericaldisp": 2,
|
||||
"mingridsampling": 3,
|
||||
}
|
||||
|
||||
# General information to configure materials
|
||||
# maxpoles: Maximum number of dispersive material poles in a model.
|
||||
@@ -330,6 +334,26 @@ class SimulationConfig:
|
||||
self._set_input_file_path()
|
||||
self._set_model_start_end()
|
||||
|
||||
def set_model_device(self, deviceID):
|
||||
"""Specify pycuda/pyopencl object for model.
|
||||
|
||||
Args:
|
||||
deviceID: int of requested deviceID of compute device.
|
||||
|
||||
Returns:
|
||||
dev: requested pycuda/pyopencl device object.
|
||||
"""
|
||||
|
||||
found = False
|
||||
for ID, dev in self.devices["devs"].items():
|
||||
if ID == deviceID:
|
||||
found = True
|
||||
return dev
|
||||
|
||||
if not found:
|
||||
logger.exception(f"Compute device with device ID {deviceID} does not exist.")
|
||||
raise ValueError
|
||||
|
||||
def _set_precision(self):
|
||||
"""Data type (precision) for electromagnetic field output.
|
||||
|
||||
@@ -389,26 +413,6 @@ class SimulationConfig:
|
||||
self.model_start = modelstart
|
||||
self.model_end = modelend
|
||||
|
||||
def get_model_device(self, deviceID):
|
||||
"""Specify pycuda/pyopencl object for model.
|
||||
|
||||
Args:
|
||||
deviceID: int of requested deviceID of compute device.
|
||||
|
||||
Returns:
|
||||
dev: requested pycuda/pyopencl device object.
|
||||
"""
|
||||
|
||||
found = False
|
||||
for ID, dev in self.devices["devs"].items():
|
||||
if ID == deviceID:
|
||||
found = True
|
||||
return dev
|
||||
|
||||
if not found:
|
||||
logger.error(f"Compute device with device ID {deviceID} does " "not exist.")
|
||||
raise ValueError
|
||||
|
||||
def get_model_config(self, model_num: Optional[int] = None) -> ModelConfig:
|
||||
"""Return ModelConfig instance for specific model.
|
||||
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -21,7 +21,6 @@ import numpy as np
|
||||
cimport numpy as np
|
||||
|
||||
np.seterr(divide='raise')
|
||||
from cython.parallel import prange
|
||||
|
||||
from gprMax.cython.yee_cell_setget_rigid cimport (
|
||||
set_rigid_E,
|
||||
@@ -686,7 +685,6 @@ cpdef void build_box(
|
||||
int yf,
|
||||
int zs,
|
||||
int zf,
|
||||
int nthreads,
|
||||
int numID,
|
||||
int numIDx,
|
||||
int numIDy,
|
||||
@@ -701,7 +699,6 @@ cpdef void build_box(
|
||||
|
||||
Args:
|
||||
xs, xf, ys, yf, zs, zf: ints for cell coordinates of entire box.
|
||||
nthreads: int for number of threads to use
|
||||
numID, numIDx, numIDy, numIDz: ints for numeric ID of material.
|
||||
averaging: bint for whether material property averaging will occur for
|
||||
the object.
|
||||
@@ -711,14 +708,14 @@ cpdef void build_box(
|
||||
cdef Py_ssize_t i, j, k
|
||||
|
||||
if averaging:
|
||||
for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for i in range(xs, xf):
|
||||
for j in range(ys, yf):
|
||||
for k in range(zs, zf):
|
||||
solid[i, j, k] = numID
|
||||
unset_rigid_E(i, j, k, rigidE)
|
||||
unset_rigid_H(i, j, k, rigidH)
|
||||
else:
|
||||
for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for i in range(xs, xf):
|
||||
for j in range(ys, yf):
|
||||
for k in range(zs, zf):
|
||||
solid[i, j, k] = numID
|
||||
@@ -731,32 +728,32 @@ cpdef void build_box(
|
||||
ID[4, i, j, k] = numIDy
|
||||
ID[5, i, j, k] = numIDz
|
||||
|
||||
for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for i in range(xs, xf):
|
||||
j = yf
|
||||
k = zf
|
||||
ID[0, i, j, k] = numIDx
|
||||
|
||||
i = xf
|
||||
for j in prange(ys, yf, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for j in range(ys, yf):
|
||||
for k in range(zf, zf + 1):
|
||||
ID[1, i, j, k] = numIDy
|
||||
|
||||
i = xf
|
||||
j = yf
|
||||
for k in prange(zs, zf, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(zs, zf):
|
||||
ID[2, i, j, k] = numIDz
|
||||
|
||||
i = xf
|
||||
for j in prange(ys, yf, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for j in range(ys, yf):
|
||||
for k in range(zs, zf):
|
||||
ID[3, i, j, k] = numIDx
|
||||
|
||||
for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for i in range(xs, xf):
|
||||
j = yf
|
||||
for k in range(zs, zf):
|
||||
ID[4, i, j, k] = numIDy
|
||||
|
||||
for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for i in range(xs, xf):
|
||||
for j in range(ys, yf):
|
||||
k = zf
|
||||
ID[5, i, j, k] = numIDz
|
||||
@@ -1316,7 +1313,6 @@ cpdef void build_voxels_from_array(
|
||||
int xs,
|
||||
int ys,
|
||||
int zs,
|
||||
int nthreads,
|
||||
int numexistmaterials,
|
||||
bint averaging,
|
||||
np.int16_t[:, :, ::1] data,
|
||||
@@ -1330,7 +1326,6 @@ cpdef void build_voxels_from_array(
|
||||
Args:
|
||||
xs, ys, zs: ints for cell coordinates of position of start of array in
|
||||
domain.
|
||||
nthreads: int for number of threads to use
|
||||
numexistmaterials: int for number of existing materials in model prior
|
||||
to building voxels.
|
||||
averaging: bint for whether material property averaging will occur for
|
||||
@@ -1364,7 +1359,7 @@ cpdef void build_voxels_from_array(
|
||||
else:
|
||||
zf = zs + data.shape[2]
|
||||
|
||||
for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for i in range(xs, xf):
|
||||
for j in range(ys, yf):
|
||||
for k in range(zs, zf):
|
||||
numID = data[i - xs, j - ys, k - zs]
|
||||
@@ -1377,7 +1372,6 @@ cpdef void build_voxels_from_array_mask(
|
||||
int xs,
|
||||
int ys,
|
||||
int zs,
|
||||
int nthreads,
|
||||
int waternumID,
|
||||
int grassnumID,
|
||||
bint averaging,
|
||||
@@ -1392,7 +1386,6 @@ cpdef void build_voxels_from_array_mask(
|
||||
|
||||
Args:
|
||||
xs, ys, zs: ints for cell coordinates of position of start of array in domain.
|
||||
nthreads: int for number of threads to use
|
||||
waternumID, grassnumID: ints for numeric ID of water and grass materials.
|
||||
averaging: bint for whether material property averaging will occur for
|
||||
the object.
|
||||
@@ -1409,7 +1402,7 @@ cpdef void build_voxels_from_array_mask(
|
||||
yf = ys + data.shape[1]
|
||||
zf = zs + data.shape[2]
|
||||
|
||||
for i in prange(xs, xf, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for i in range(xs, xf):
|
||||
for j in range(ys, yf):
|
||||
for k in range(zs, zf):
|
||||
if mask[i - xs, j - ys, k - zs] == 1:
|
||||
|
||||
892
gprMax/cython/plane_wave.pyx
普通文件
892
gprMax/cython/plane_wave.pyx
普通文件
@@ -0,0 +1,892 @@
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, John Hartley, and
|
||||
# Adittya Pal
|
||||
#
|
||||
# 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
|
||||
|
||||
cimport cython
|
||||
from libc.math cimport M_PI, abs, atan2, ceil, cos, exp, floor, pow, round, sin, sqrt, tan
|
||||
from libc.stdio cimport FILE, fclose, fopen, fwrite
|
||||
from libc.string cimport strcmp
|
||||
|
||||
from cython.parallel import prange
|
||||
|
||||
from gprMax.config cimport float_or_double
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cpdef double[:, ::1] getProjections(
|
||||
double psi,
|
||||
int[:] m
|
||||
):
|
||||
"""Gets the projection vectors to source magnetic fields of plane wave.
|
||||
|
||||
Args:
|
||||
psi: float for angle describing polatan value of required phi angle
|
||||
(which would be approximated to a rational number).
|
||||
m: int array to store integer mappings, m_x, m_y, m_z which determine
|
||||
the rational angles, for assignment of the correct element to 3D
|
||||
FDTD grid from 1D representation, last element stores
|
||||
max(m_x, m_y, m_z).
|
||||
|
||||
Returns:
|
||||
projections: float array to store projections for sourcing magnetic
|
||||
field and the sourcing vector.
|
||||
"""
|
||||
|
||||
cdef double phi, theta, cos_phi, sin_phi, cos_psi, sin_psi, cos_theta, sin_theta
|
||||
cdef double[:, ::1] projections = np.zeros((2, 3), order='C')
|
||||
|
||||
if m[0] == 0:
|
||||
phi = M_PI / 2
|
||||
else:
|
||||
phi = atan2(m[1], m[0])
|
||||
|
||||
if m[2] == 0:
|
||||
theta = M_PI / 2
|
||||
else:
|
||||
theta = atan2(sqrt(m[0] * m[0] + m[1] * m[1]), m[2])
|
||||
|
||||
cos_phi = cos(phi)
|
||||
sin_phi = sin(phi)
|
||||
cos_theta = cos(theta)
|
||||
sin_theta = sin(theta)
|
||||
cos_psi = cos(psi)
|
||||
sin_psi = sin(psi)
|
||||
|
||||
# Magnetic field projection vector
|
||||
projections[0, 0] = sin_psi * sin_phi + cos_psi * cos_theta * cos_phi
|
||||
projections[0, 1] = -sin_psi * cos_phi + cos_psi * cos_theta * sin_phi
|
||||
projections[0, 2] = -cos_psi * sin_theta
|
||||
|
||||
# Direction cosines
|
||||
projections[1, 0] = sin_theta * cos_phi
|
||||
projections[1, 1] = sin_theta * sin_phi
|
||||
projections[1, 2] = cos_theta
|
||||
|
||||
return projections
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cdef int[:] getPhi(
|
||||
int[:, :] integers,
|
||||
double required_ratio,
|
||||
double tolerance
|
||||
):
|
||||
"""Gets rational angle approximation to phi within the requested tolerance
|
||||
level using Farey Fractions to determine a rational number closest to
|
||||
the real number.
|
||||
|
||||
Args:
|
||||
integers: int array to determine the value of m_x and m_y.
|
||||
required_ratio: float of tan value of the required phi angle
|
||||
(which would be approximated to a rational number).
|
||||
tolerance: float for acceptable deviation in the tan value of the
|
||||
rational angle from phi.
|
||||
|
||||
Returns:
|
||||
integers: int array of sequence of the two integers [m_x, m_y].
|
||||
"""
|
||||
|
||||
if(abs(integers[2, 0]/<double>integers[2, 1]-required_ratio)<tolerance):
|
||||
return integers[2, :]
|
||||
while True:
|
||||
integers[1, 0] = integers[0, 0]+integers[2, 0]
|
||||
integers[1, 1] = integers[0, 1]+integers[2, 1]
|
||||
ratio = integers[1, 0]/<double>integers[1, 1]
|
||||
if (abs(ratio - required_ratio) <= tolerance):
|
||||
return integers[1, :]
|
||||
elif(ratio >= required_ratio):
|
||||
integers[2, :] = integers[1, :]
|
||||
else:
|
||||
integers[0, :] = integers[1, :]
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cdef inline double getTanValue(
|
||||
int[:] integers,
|
||||
double[:] dr
|
||||
):
|
||||
"""Returns tan value of the angle approximated to theta given three integers.
|
||||
|
||||
Args:
|
||||
integers: int array of three integers for the rational angle
|
||||
approximation.
|
||||
dr: double array containing the separation between grid points along
|
||||
the three axes [dx, dy, dz].
|
||||
|
||||
Returns:
|
||||
_tanValue: double of tan value of the rational angle corresponding to
|
||||
integers m_x, m_y, m_z.
|
||||
"""
|
||||
|
||||
if(integers[2]==0): #if rational angle==90 degrees
|
||||
return 99999.0 #return a large number to avoid division by zero error
|
||||
else:
|
||||
return sqrt((integers[0]/dr[0])*(integers[0]/dr[0]) + (integers[1]/dr[1])*(integers[1]/dr[1])
|
||||
)/(integers[2]/dr[2])
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cdef int[:, :] get_mZ(
|
||||
int m_x,
|
||||
int m_y,
|
||||
double theta,
|
||||
double[:] Delta_r
|
||||
):
|
||||
"""Gets arrays to perform a binary search to determine a rational number,
|
||||
m_z, closest to real number, m_z, to get desired tan Theta value.
|
||||
|
||||
Args:
|
||||
m_x and m_y: ints approximating rational angle to tan value of phi.
|
||||
theta: float of polar angle of incident plane wave (radians) to be
|
||||
approximated to a rational angle.
|
||||
Delta_r: float array of projections of propagation vector along
|
||||
different coordinate axes.
|
||||
|
||||
Returns:
|
||||
_integers: int array of 2D sequence of three integers [m_x, m_y, m_z]
|
||||
to perform a binary search to determine value of m_z within
|
||||
given limits.
|
||||
"""
|
||||
|
||||
cdef double m_z = 0
|
||||
m_z = sqrt((m_x/Delta_r[0])*(m_x/Delta_r[0]) + (m_y/Delta_r[1])*(m_y/Delta_r[1]))/(tan(theta)/Delta_r[2]) #get an estimate of the m_z value
|
||||
return np.array([[m_x, m_y, floor(m_z)],
|
||||
[m_x, m_y, round(m_z)],
|
||||
[m_x, m_y, ceil(m_z)]], dtype=np.int32, order='C') #set up the integer array to search for an appropriate m_z
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cdef int[:] getTheta(
|
||||
int m_x,
|
||||
int m_y,
|
||||
double theta,
|
||||
double Delta_theta,
|
||||
double[:] Delta_r
|
||||
):
|
||||
"""Gets rational angle approximation to theta within requested tolerance
|
||||
level using Binary Search to determine a rational number closest to
|
||||
real number.
|
||||
|
||||
Args:
|
||||
m_x and m_y: ints approximating rational angle to tan value of phi.
|
||||
theta: float of polar angle of incident plane wave (radians) to be
|
||||
approximated to a rational angle.
|
||||
Delta_theta: float of permissible error in rational angle approximation
|
||||
to theta (radians).
|
||||
Delta_r: float array of projections of propagation vector along
|
||||
different coordinate axes.
|
||||
|
||||
Returns:
|
||||
integers: int array of sequence of three integers [m_x, m_y, m_z].
|
||||
"""
|
||||
|
||||
cdef Py_ssize_t i, j = 0
|
||||
cdef double tan_theta = 0.0
|
||||
cdef int[:, :] integers = get_mZ(m_x, m_y, theta, Delta_r) #set up the integer array to search for an appropriate m_z
|
||||
while True: #if tan value of m_z greater than permitted tolerance
|
||||
tan_theta = getTanValue(integers[1, :], Delta_r)
|
||||
if(abs(tan_theta - tan(theta)) <= Delta_theta / (cos(theta) * cos(theta))):
|
||||
break
|
||||
for i in range(3):
|
||||
for j in range(3):
|
||||
integers[i, j] = 2*integers[i, j] #expand the serach space by multiplying 2 to both the numerator and the denominator
|
||||
while(integers[0, 2]<integers[2, 2]-1): #while there are integers to search for in the denominator
|
||||
integers[1, 2] = (integers[0, 2]+integers[2, 2])/2 #get the mid integer between the upper and lower limits of denominator
|
||||
tan_theta = getTanValue(integers[1, :], Delta_r)
|
||||
if(tan_theta < tan(theta)): #if m_z results in a smaller tan value, make the denominator smaller
|
||||
integers[2, 2] = integers[1, 2] #decrease m_z, serach in the lower half of the sample space
|
||||
elif(tan_theta > tan(theta)): #if m_z results in a larger tan value, make the denominator larger
|
||||
integers[0, 2] = integers[1, 2] #increase m_z, serach in the upper half of the sample space
|
||||
|
||||
return integers[1, :]
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cpdef int[:, ::1] getIntegerForAngles(
|
||||
double phi,
|
||||
double Delta_phi,
|
||||
double theta,
|
||||
double Delta_theta,
|
||||
double[:] Delta_r
|
||||
):
|
||||
"""Gets [m_x, m_y, m_z] to determine rational angles given phi and theta
|
||||
along with the permissible tolerances.
|
||||
|
||||
Args:
|
||||
phi: float of azimuthal angle of incident plane wave (degrees) to be
|
||||
approximated to a rational angle.
|
||||
Delta_phi: float of permissible error in rational angle approximation
|
||||
to phi (degrees).
|
||||
theta: float of polar angle of incident plane wave (degrees) to be
|
||||
approximated to a rational angle.
|
||||
Delta_theta: float of permissible error in rational angle approximation
|
||||
to theta (degrees).
|
||||
Delta_r: float of projections of propagation vector along different
|
||||
coordinate axes.
|
||||
|
||||
Returns:
|
||||
quadrants[0, :]: int array specifying direction of propagation of plane
|
||||
wave along the three coordinate axes.
|
||||
quadrants[1, :]: int array of three integers [m_x, m_y, m_z].
|
||||
"""
|
||||
|
||||
cdef double required_ratio_phi, tolerance_phi = 0.0
|
||||
cdef int m_x, m_y, m_z = 0
|
||||
cdef int[:, ::1] quadrants = np.ones((2, 3), dtype=np.int32)
|
||||
if(theta>=90):
|
||||
quadrants[0, 2] = -1
|
||||
theta = 180-theta
|
||||
if(phi>=90 and phi<180):
|
||||
quadrants[0, 0] = -1
|
||||
phi = 180-phi
|
||||
elif(phi>=180 and phi<270):
|
||||
quadrants[0, 0] = -1
|
||||
quadrants[0, 1] = -1
|
||||
phi = phi-180
|
||||
elif(phi>=270 and phi<360):
|
||||
quadrants[0, 1] = -1
|
||||
phi = 360-phi
|
||||
|
||||
|
||||
if(0 <= phi < 90): #handle the case of phi=90 degrees separately
|
||||
required_ratio_phi = tan(M_PI/180*phi) * Delta_r[1] / Delta_r[0] #to avoid division by zero exception
|
||||
tolerance_phi = M_PI/180*Delta_phi / (cos(M_PI/180*phi)*cos(M_PI/180*phi)) * Delta_r[1] / Delta_r[0] #get the persissible error in tan phi
|
||||
m_y, m_x = getPhi(np.array([[floor(required_ratio_phi), 1],
|
||||
[required_ratio_phi, 1],
|
||||
[ceil(required_ratio_phi), 1]], dtype=np.int32, order='C')
|
||||
, required_ratio_phi, tolerance_phi) #get the integers [m_x, m_y] for rational angle approximation to phi
|
||||
else:
|
||||
m_x = 0
|
||||
m_y = 1
|
||||
if(theta < 90):
|
||||
m_x, m_y, m_z = getTheta(m_x, m_y, M_PI/180*theta, M_PI/180*Delta_theta, Delta_r) #get the integers [m_x, m_y, m_z] for rational angle approximation to theta
|
||||
else: #handle the case of theta=90 degrees separately
|
||||
m_z = 0 #to avoid division by zero exception
|
||||
quadrants[1, 0] = m_x
|
||||
quadrants[1, 1] = m_y
|
||||
quadrants[1, 2] = m_z
|
||||
return quadrants
|
||||
|
||||
|
||||
@cython.wraparound(False)
|
||||
@cython.boundscheck(False)
|
||||
cdef void applyTFSFMagnetic(
|
||||
int nthreads,
|
||||
float_or_double[:, :, ::1] Hx,
|
||||
float_or_double[:, :, ::1] Hy,
|
||||
float_or_double[:, :, ::1] Hz,
|
||||
float_or_double[:, ::1] E_fields,
|
||||
float_or_double[:] updatecoeffsH,
|
||||
int[:] m,
|
||||
int[:] corners
|
||||
):
|
||||
"""Implements total field-scattered field formulation for magnetic field on
|
||||
the edge of the TF/SF region of the TFSF Box.
|
||||
|
||||
Args:
|
||||
nthreads: int of number of threads to parallelize for loops.
|
||||
Hx, Hy, Hz: double array to store magnetic fields for grid cells over
|
||||
the TFSF box at particular indices.
|
||||
E_fields: double array to store electric fields of 1D representation of
|
||||
plane wave in a direction along which the wave propagates.
|
||||
updatecoeffsH: float of coefficients of fields in TFSF assignment
|
||||
equation for the magnetic field.
|
||||
m: int array of integer mappings, m_x, m_y, m_z which determine rational
|
||||
angles for assignment of correct element to 3D FDTD grid from 1D
|
||||
representation, last element stores max(m_x, m_y, m_z).
|
||||
corners: int array of coordinates of corners of TF/SF field boundaries.
|
||||
"""
|
||||
|
||||
cdef Py_ssize_t i, j, k = 0
|
||||
|
||||
# Precompute index values
|
||||
cdef int index = 0
|
||||
cdef int m_x = m[0]
|
||||
cdef int m_y = m[1]
|
||||
cdef int m_z = m[2]
|
||||
|
||||
cdef int x_start = corners[0]
|
||||
cdef int y_start = corners[1]
|
||||
cdef int z_start = corners[2]
|
||||
cdef int x_stop = corners[3]
|
||||
cdef int y_stop = corners[4]
|
||||
cdef int z_stop = corners[5]
|
||||
|
||||
cdef float_or_double[:] E_x = E_fields[0, :]
|
||||
cdef float_or_double[:] E_y = E_fields[1, :]
|
||||
cdef float_or_double[:] E_z = E_fields[2, :]
|
||||
|
||||
cdef float_or_double coef_H_xy = updatecoeffsH[2]
|
||||
cdef float_or_double coef_H_xz = updatecoeffsH[3]
|
||||
cdef float_or_double coef_H_yz = updatecoeffsH[3]
|
||||
cdef float_or_double coef_H_yx = updatecoeffsH[1]
|
||||
cdef float_or_double coef_H_zx = updatecoeffsH[1]
|
||||
cdef float_or_double coef_H_zy = updatecoeffsH[2]
|
||||
|
||||
#**** constant x faces -- scattered-field nodes ****
|
||||
i = x_start
|
||||
for j in prange(y_start, y_stop+1, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop):
|
||||
#correct Hy at firstX-1/2 by subtracting Ez_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Hy[i-1, j, k] -= coef_H_yx * E_z[index]
|
||||
|
||||
for j in prange(y_start, y_stop, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop+1):
|
||||
#correct Hz at firstX-1/2 by adding Ey_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Hz[i-1, j, k] += coef_H_zx * E_y[index]
|
||||
|
||||
i = x_stop
|
||||
for j in prange(y_start, y_stop+1, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop):
|
||||
#correct Hy at lastX+1/2 by adding Ez_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Hy[i, j, k] += coef_H_yx * E_z[index]
|
||||
|
||||
for j in prange(y_start, y_stop, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop+1):
|
||||
#correct Hz at lastX+1/2 by subtractinging Ey_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Hz[i, j, k] -= coef_H_zx * E_y[index]
|
||||
|
||||
#**** constant y faces -- scattered-field nodes ****
|
||||
j = y_start
|
||||
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop):
|
||||
#correct Hx at firstY-1/2 by adding Ez_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Hx[i, j-1, k] += coef_H_xy * E_z[index]
|
||||
|
||||
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop+1):
|
||||
#correct Hz at firstY-1/2 by subtracting Ex_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Hz[i, j-1, k] -= coef_H_zy * E_x[index]
|
||||
|
||||
j = y_stop
|
||||
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop):
|
||||
#correct Hx at lastY+1/2 by subtracting Ez_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Hx[i, j, k] -= coef_H_xy * E_z[index]
|
||||
|
||||
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop+1):
|
||||
#correct Hz at lastY-1/2 by adding Ex_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Hz[i, j, k] += coef_H_zy * E_x[index]
|
||||
|
||||
#**** constant z faces -- scattered-field nodes ****
|
||||
k = z_start
|
||||
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for j in range(y_start, y_stop+1):
|
||||
#correct Hy at firstZ-1/2 by adding Ex_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Hy[i, j, k-1] += coef_H_yz * E_x[index]
|
||||
|
||||
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for j in range(y_start, y_stop):
|
||||
#correct Hx at firstZ-1/2 by subtracting Ey_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Hx[i, j, k-1] -= coef_H_xz * E_y[index]
|
||||
|
||||
k = z_stop
|
||||
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for j in range(y_start, y_stop+1):
|
||||
#correct Hy at firstZ-1/2 by subtracting Ex_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Hy[i, j, k] -= coef_H_yz * E_x[index]
|
||||
|
||||
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for j in range(y_start, y_stop):
|
||||
#correct Hx at lastZ+1/2 by adding Ey_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Hx[i, j, k] += coef_H_xz * E_y[index]
|
||||
|
||||
|
||||
cdef void applyTFSFElectric(
|
||||
int nthreads,
|
||||
float_or_double[:, :, ::1] Ex,
|
||||
float_or_double[:, :, ::1] Ey,
|
||||
float_or_double[:, :, ::1] Ez,
|
||||
float_or_double[:, ::1] H_fields,
|
||||
float_or_double[:] updatecoeffsE,
|
||||
int[:] m,
|
||||
int[:] corners
|
||||
):
|
||||
|
||||
"""Implements total field-scattered field formulation for electric field on
|
||||
edge of the TF/SF region of the TFSF Box.
|
||||
|
||||
Args:
|
||||
nthreads: int for number of threads to parallelize the for loops.
|
||||
Ex, Ey, Ez: double array for magnetic fields for grid cells over TFSF
|
||||
box at particular indices.
|
||||
H_fields: double array to store electric fields of 1D representation of
|
||||
plane wave in direction along which wave propagates.
|
||||
updatecoeffsE: float of coefficients of fields in TFSF assignment
|
||||
equation for magnetic field.
|
||||
m: int array of integer mappings, m_x, m_y, m_z which determine rational
|
||||
angles for assignment of correct element to 3D FDTD grid from 1D
|
||||
representation, last element stores max(m_x, m_y, m_z).
|
||||
corners: int array for coordinates of corners of TF/SF field boundaries.
|
||||
"""
|
||||
|
||||
cdef Py_ssize_t i, j, k = 0
|
||||
|
||||
# Precompute index values
|
||||
cdef int index = 0
|
||||
cdef int m_x = m[0]
|
||||
cdef int m_y = m[1]
|
||||
cdef int m_z = m[2]
|
||||
|
||||
cdef int x_start = corners[0]
|
||||
cdef int y_start = corners[1]
|
||||
cdef int z_start = corners[2]
|
||||
cdef int x_stop = corners[3]
|
||||
cdef int y_stop = corners[4]
|
||||
cdef int z_stop = corners[5]
|
||||
|
||||
cdef float_or_double[:] H_x = H_fields[0, :]
|
||||
cdef float_or_double[:] H_y = H_fields[1, :]
|
||||
cdef float_or_double[:] H_z = H_fields[2, :]
|
||||
|
||||
cdef float_or_double coef_E_xy = updatecoeffsE[2]
|
||||
cdef float_or_double coef_E_xz = updatecoeffsE[3]
|
||||
cdef float_or_double coef_E_yz = updatecoeffsE[3]
|
||||
cdef float_or_double coef_E_yx = updatecoeffsE[1]
|
||||
cdef float_or_double coef_E_zx = updatecoeffsE[1]
|
||||
cdef float_or_double coef_E_zy = updatecoeffsE[2]
|
||||
|
||||
#**** constant x faces -- total-field nodes ****/
|
||||
i = x_start
|
||||
for j in prange(y_start, y_stop+1, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop):
|
||||
#correct Ez at firstX face by subtracting Hy_inc
|
||||
index = m_x * (i-1) + m_y * j + m_z * k
|
||||
Ez[i, j, k] -= coef_E_zx * H_y[index]
|
||||
|
||||
for j in prange(y_start, y_stop, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop+1):
|
||||
#correct Ey at firstX face by adding Hz_inc
|
||||
index = m_x * (i-1) + m_y * j + m_z * k
|
||||
Ey[i, j, k] += coef_E_yx * H_z[index]
|
||||
|
||||
i = x_stop
|
||||
for j in prange(y_start, y_stop+1, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop):
|
||||
#correct Ez at lastX face by adding Hy_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Ez[i, j, k] += coef_E_zx * H_y[index]
|
||||
|
||||
i = x_stop
|
||||
for j in prange(y_start, y_stop, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop+1):
|
||||
#correct Ey at lastX face by subtracting Hz_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Ey[i, j, k] -= coef_E_yx * H_z[index]
|
||||
|
||||
#**** constant y faces -- total-field nodes ****/
|
||||
j = y_start
|
||||
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop):
|
||||
#correct Ez at firstY face by adding Hx_inc
|
||||
index = m_x * i + m_y * (j-1) + m_z * k
|
||||
Ez[i, j, k] += coef_E_zy * H_x[index]
|
||||
|
||||
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop+1):
|
||||
#correct Ex at firstY face by subtracting Hz_inc
|
||||
index = m_x * i + m_y * (j-1) + m_z * k
|
||||
Ex[i, j, k] -= coef_E_xy * H_z[index]
|
||||
|
||||
j = y_stop
|
||||
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop):
|
||||
#correct Ez at lastY face by subtracting Hx_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Ez[i, j, k] -= coef_E_zy * H_x[index]
|
||||
|
||||
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for k in range(z_start, z_stop+1):
|
||||
#correct Ex at lastY face by adding Hz_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Ex[i, j, k] += coef_E_xy * H_z[index]
|
||||
|
||||
#**** constant z faces -- total-field nodes ****/
|
||||
k = z_start
|
||||
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for j in range(y_start, y_stop):
|
||||
#correct Ey at firstZ face by subtracting Hx_inc
|
||||
index = m_x * i + m_y * j + m_z * (k-1)
|
||||
Ey[i, j, k] -= coef_E_yz * H_x[index]
|
||||
|
||||
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for j in range(y_start, y_stop+1):
|
||||
#correct Ex at firstZ face by adding Hy_inc
|
||||
index = m_x * i + m_y * j + m_z * (k-1)
|
||||
Ex[i, j, k] += coef_E_xz * H_y[index]
|
||||
|
||||
k = z_stop
|
||||
for i in prange(x_start, x_stop+1, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for j in range(y_start, y_stop):
|
||||
#correct Ey at lastZ face by adding Hx_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Ey[i, j, k] += coef_E_yz * H_x[index]
|
||||
|
||||
for i in prange(x_start, x_stop, nogil=True, schedule='static', num_threads=nthreads):
|
||||
for j in range(y_start, y_stop+1):
|
||||
#correct Ex at lastZ face by subtracting Hy_inc
|
||||
index = m_x * i + m_y * j + m_z * k
|
||||
Ex[i, j, k] -= coef_E_xz * H_y[index]
|
||||
|
||||
|
||||
cdef void initializeMagneticFields(
|
||||
int[:] m,
|
||||
float_or_double[:, ::1] H_fields,
|
||||
double[:] projections,
|
||||
float_or_double[:, ::1] waveformvalues_wholedt,
|
||||
bint precompute,
|
||||
int iteration,
|
||||
double dt,
|
||||
double ds,
|
||||
double c,
|
||||
double start,
|
||||
double stop,
|
||||
double freq,
|
||||
char* wavetype
|
||||
):
|
||||
"""Initialises first few grid points of source waveform.
|
||||
|
||||
Args:
|
||||
m: int array of integer mappings, m_x, m_y, m_z which determine rational
|
||||
angles for assignment of correct element to 3D FDTD grid from 1D
|
||||
representation, last element stores max(m_x, m_y, m_z).
|
||||
H_fields: double array for electric fields of 1D representation of plane
|
||||
wave in a direction along which the wave propagates.
|
||||
projections: float array for projections of magnetic fields along
|
||||
different cartesian axes.
|
||||
waveformvalues_wholedt: double array stores precomputed waveforms at
|
||||
each timestep to initialize magnetic fields.
|
||||
precompute: boolean to store whether fields values have been precomputed
|
||||
or should be computed on the fly.
|
||||
iterations: int stores number of iterations in the simulation.
|
||||
dt: float of timestep for the simulation.
|
||||
ds: float of projection vector for sourcing the plane wave.
|
||||
c: float of speed of light in the medium.
|
||||
start: float of start time at which source is placed in the TFSF grid.
|
||||
stop: float of stop time at which source is removed from TFSF grid.
|
||||
freq: float of frequency of introduced wave which determines grid points
|
||||
per wavelength for wave source.
|
||||
wavetype: string of type of waveform whose magnitude should be returned.
|
||||
"""
|
||||
|
||||
cdef Py_ssize_t r = 0
|
||||
cdef double time_x, time_y, time_z = 0.0
|
||||
if (precompute == True):
|
||||
for r in range(m[3]): #loop to assign the source values of magnetic field to the first few gridpoints
|
||||
H_fields[0, r] = projections[0] * waveformvalues_wholedt[0, r]
|
||||
H_fields[1, r] = projections[1] * waveformvalues_wholedt[1, r]
|
||||
H_fields[2, r] = projections[2] * waveformvalues_wholedt[2, r]
|
||||
else:
|
||||
for r in range(m[3]): #loop to assign the source values of magnetic field to the first few gridpoints
|
||||
time_x = dt * iteration - (r+ (m[1]+m[2])*0.5) * ds/c
|
||||
time_y = dt * iteration - (r+ (m[2]+m[0])*0.5) * ds/c
|
||||
time_z = dt * iteration - (r+ (m[0]+m[1])*0.5) * ds/c
|
||||
if (dt * iteration >= start and dt * iteration <= stop):
|
||||
# Set the time of the waveform evaluation to account for any delay in the start
|
||||
H_fields[0, r] = projections[0] * getSource(time_x-start, freq, wavetype, dt)
|
||||
H_fields[1, r] = projections[1] * getSource(time_y-start, freq, wavetype, dt)
|
||||
H_fields[2, r] = projections[2] * getSource(time_z-start, freq, wavetype, dt)
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cdef void updateMagneticFields(
|
||||
int n,
|
||||
float_or_double[:, ::1] H_fields,
|
||||
float_or_double[:, ::1] E_fields,
|
||||
float_or_double[:] updatecoeffsH,
|
||||
int[:] m
|
||||
):
|
||||
"""Updates magnetic fields for next time step using Equation 8 of
|
||||
DOI: 10.1109/LAWP.2009.2016851.
|
||||
|
||||
Args:
|
||||
n: int for spatial length of DPW array to update each length grid
|
||||
cell.
|
||||
H_fields: double array of magnetic fields of DPW until temporal
|
||||
index time.
|
||||
E_fields: double array of electric fields of DPW until temporal
|
||||
index time.
|
||||
updatecoeffsH: double array of coefficients of fields in update
|
||||
equation for magnetic field.
|
||||
m: int array of integer mappings, m_x, m_y, m_z which determine
|
||||
rational angles for assignment of correct element to 3D FDTD
|
||||
grid from 1D representation, last element stores
|
||||
max(m_x, m_y, m_z).
|
||||
"""
|
||||
|
||||
cdef Py_ssize_t j = 0
|
||||
|
||||
cdef float_or_double[:] E_x = E_fields[0, :]
|
||||
cdef float_or_double[:] E_y = E_fields[1, :]
|
||||
cdef float_or_double[:] E_z = E_fields[2, :]
|
||||
cdef float_or_double[:] H_x = H_fields[0, :]
|
||||
cdef float_or_double[:] H_y = H_fields[1, :]
|
||||
cdef float_or_double[:] H_z = H_fields[2, :]
|
||||
|
||||
cdef float_or_double coef_H_xt = updatecoeffsH[0]
|
||||
cdef float_or_double coef_H_xy = updatecoeffsH[2]
|
||||
cdef float_or_double coef_H_xz = updatecoeffsH[3]
|
||||
|
||||
cdef float_or_double coef_H_yt = updatecoeffsH[0]
|
||||
cdef float_or_double coef_H_yz = updatecoeffsH[3]
|
||||
cdef float_or_double coef_H_yx = updatecoeffsH[1]
|
||||
|
||||
cdef float_or_double coef_H_zt = updatecoeffsH[0]
|
||||
cdef float_or_double coef_H_zx = updatecoeffsH[1]
|
||||
cdef float_or_double coef_H_zy = updatecoeffsH[2]
|
||||
|
||||
cdef int m_x = m[0]
|
||||
cdef int m_y = m[1]
|
||||
cdef int m_z = m[2]
|
||||
|
||||
for j in range(m[3], n-m[3]): #loop to update the magnetic field at each spatial index
|
||||
H_x[j] = coef_H_xt * H_x[j] + coef_H_xz * ( E_y[j+m_z] - E_y[j] ) - coef_H_xy * ( E_z[j+m_y] - E_z[j] ) #equation 8 of Tan, Potter paper
|
||||
H_y[j] = coef_H_yt * H_y[j] + coef_H_yx * ( E_z[j+m_x] - E_z[j] ) - coef_H_yz * ( E_x[j+m_z] - E_x[j] ) #equation 8 of Tan, Potter paper
|
||||
H_z[j] = coef_H_zt * H_z[j] + coef_H_zy * ( E_x[j+m_y] - E_x[j] ) - coef_H_zx * ( E_y[j+m_x] - E_y[j] ) #equation 8 of Tan, Potter paper
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cdef void updateElectricFields(
|
||||
int n,
|
||||
float_or_double[:, ::1] H_fields,
|
||||
float_or_double[:, ::1] E_fields,
|
||||
float_or_double[:] updatecoeffsE,
|
||||
int[:] m
|
||||
):
|
||||
"""Updates electric fields for next time step using Equation 9 of
|
||||
DOI: 10.1109/LAWP.2009.2016851.
|
||||
|
||||
Args:
|
||||
n: int for spatial length of DPW array to update each length grid
|
||||
cell.
|
||||
H_fields: double array of magnetic fields of DPW until temporal
|
||||
index time.
|
||||
E_fields: double array of electric fields of DPW until temporal
|
||||
index time.
|
||||
updatecoeffsE: double array of coefficients of fields in update
|
||||
equation for electric field.
|
||||
m: int array of integer mappings, m_x, m_y, m_z which determine
|
||||
rational angles for assignment of correct element to 3D FDTD
|
||||
grid from 1D representation, last element stores
|
||||
max(m_x, m_y, m_z).
|
||||
"""
|
||||
|
||||
cdef Py_ssize_t j = 0
|
||||
|
||||
cdef float_or_double[:] E_x = E_fields[0, :]
|
||||
cdef float_or_double[:] E_y = E_fields[1, :]
|
||||
cdef float_or_double[:] E_z = E_fields[2, :]
|
||||
cdef float_or_double[:] H_x = H_fields[0, :]
|
||||
cdef float_or_double[:] H_y = H_fields[1, :]
|
||||
cdef float_or_double[:] H_z = H_fields[2, :]
|
||||
|
||||
cdef float_or_double coef_E_xt = updatecoeffsE[0]
|
||||
cdef float_or_double coef_E_xy = updatecoeffsE[2]
|
||||
cdef float_or_double coef_E_xz = updatecoeffsE[3]
|
||||
|
||||
cdef float_or_double coef_E_yt = updatecoeffsE[0]
|
||||
cdef float_or_double coef_E_yz = updatecoeffsE[3]
|
||||
cdef float_or_double coef_E_yx = updatecoeffsE[1]
|
||||
|
||||
cdef float_or_double coef_E_zt = updatecoeffsE[0]
|
||||
cdef float_or_double coef_E_zx = updatecoeffsE[1]
|
||||
cdef float_or_double coef_E_zy = updatecoeffsE[2]
|
||||
|
||||
cdef int m_x = m[0]
|
||||
cdef int m_y = m[1]
|
||||
cdef int m_z = m[2]
|
||||
|
||||
for j in range(m[3], n-m[3]): #loop to update the electric field at each spatial index
|
||||
E_x[j] = coef_E_xt * E_x[j] + coef_E_xz * ( H_z[j] - H_z[j-m_y] ) - coef_E_xy * ( H_y[j] - H_y[j-m_z] ) #equation 9 of Tan, Potter paper
|
||||
E_y[j] = coef_E_yt * E_y[j] + coef_E_yx * ( H_x[j] - H_x[j-m_z] ) - coef_E_yz * ( H_z[j] - H_z[j-m_x] ) #equation 9 of Tan, Potter paper
|
||||
E_z[j] = coef_E_zt * E_z[j] + coef_E_zy * ( H_y[j] - H_y[j-m_x] ) - coef_E_zx * ( H_x[j] - H_x[j-m_y] ) #equation 9 of Tan, Potter paper
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cpdef double getSource(
|
||||
double time,
|
||||
double freq,
|
||||
char* wavetype,
|
||||
double dt
|
||||
):
|
||||
"""Gets magnitude of source field in direction perpendicular to propagation
|
||||
of plane wave.
|
||||
|
||||
Args:
|
||||
time: float of time at which magnitude of source is calculated.
|
||||
freq: float of frequency of introduced wave which determines grid points
|
||||
per wavelength for wave source.
|
||||
wavetype: string of type of waveform whose magnitude should be returned.
|
||||
dt: double of time upto which wave should exist in a impulse delta pulse.
|
||||
|
||||
Returns:
|
||||
sourceMagnitude: double of magnitude of source for requested indices at
|
||||
current time.
|
||||
"""
|
||||
|
||||
# Waveforms
|
||||
if (strcmp(wavetype, "gaussian") == 0):
|
||||
return exp(-2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0)))
|
||||
|
||||
elif (strcmp(wavetype, "gaussiandot") == 0 or strcmp(wavetype, "gaussianprime") == 0):
|
||||
return -4.0 * M_PI * M_PI * freq * (time * freq - 1.0
|
||||
) * exp(-2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0)))
|
||||
|
||||
elif (strcmp(wavetype, "gaussiandotnorm") == 0):
|
||||
return -2.0 * M_PI * (time * freq - 1.0
|
||||
) * exp(-2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0))) * exp(0.5)
|
||||
|
||||
elif (strcmp(wavetype, "gaussiandotdot") == 0 or strcmp(wavetype, "gaussiandoubleprime") == 0):
|
||||
return (2.0 * M_PI * freq) * (2.0 * M_PI * freq) * (2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0)) - 1.0
|
||||
) * exp(-2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0)))
|
||||
|
||||
elif (strcmp(wavetype, "gaussiandotdotnorm") == 0):
|
||||
return (2.0 * (M_PI *(time * freq - 1.0)) * (M_PI * (time * freq - 1.0)) - 1.0
|
||||
) * exp(-2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0)))
|
||||
|
||||
elif (strcmp(wavetype, "ricker") == 0):
|
||||
return (1.0 - 2.0 * (M_PI *(time * freq - 1.0)) * (M_PI * (time * freq - 1.0))
|
||||
) * exp(-2.0 * (M_PI * (time * freq - 1.0)) * (M_PI * (time * freq - 1.0))) # define a Ricker wave source
|
||||
|
||||
elif (strcmp(wavetype, "sine") == 0):
|
||||
if (time * freq <= 1):
|
||||
return sin(2.0 * M_PI * freq * time)
|
||||
else:
|
||||
return 0.0
|
||||
|
||||
elif (strcmp(wavetype, "contsine") == 0):
|
||||
return min(0.25 * time* freq, 1) * sin(2 * M_PI * time* freq)
|
||||
|
||||
elif (strcmp(wavetype, "impulse") == 0):
|
||||
if (time < dt): # time < dt condition required to do impulsive magnetic dipole
|
||||
return 1.0
|
||||
else:
|
||||
return 0.0
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cpdef void calculate1DWaveformValues(
|
||||
float_or_double[:, :, ::1] waveformvalues_wholedt,
|
||||
int iterations,
|
||||
int[:] m,
|
||||
double dt,
|
||||
double ds,
|
||||
double c,
|
||||
double start,
|
||||
double stop,
|
||||
double freq,
|
||||
char* wavetype
|
||||
):
|
||||
"""Precomputes source waveform values so that initialization is faster,
|
||||
if requested.
|
||||
|
||||
Args:
|
||||
waveformvalues_wholedt: double array of precomputed waveforms at each
|
||||
timestep to initialize magnetic fields.
|
||||
iterations: int of number of iterations in simulation.
|
||||
m: int array of integer mappings, m_x, m_y, m_z which determine rational
|
||||
angles for assignment of correct element to 3D FDTD grid from 1D
|
||||
representation, last element stores max(m_x, m_y, m_z).
|
||||
dt: float of timestep for the simulation.
|
||||
ds: float of projection vector for sourcing the plane wave.
|
||||
c: float of speed of light in the medium.
|
||||
start: float of start time at which source is placed in the TFSF grid.
|
||||
stop: float of stop time at which source is removed from TFSF grid.
|
||||
freq: float of frequency of introduced wave which determines grid points
|
||||
per wavelength for wave source.
|
||||
wavetype: string of type of waveform whose magnitude should be returned.
|
||||
"""
|
||||
|
||||
cdef double time_x, time_y, time_z = 0.0
|
||||
cdef Py_ssize_t iteration, r = 0
|
||||
|
||||
for iteration in range(iterations):
|
||||
for r in range(m[3]):
|
||||
time_x = dt * iteration - (r+ (m[1]+m[2])*0.5) * ds/c
|
||||
time_y = dt * iteration - (r+ (m[2]+m[0])*0.5) * ds/c
|
||||
time_z = dt * iteration - (r+ (m[0]+m[1])*0.5) * ds/c
|
||||
if (dt * iteration >= start and dt * iteration <= stop):
|
||||
# Set the time of the waveform evaluation to account for any delay in the start
|
||||
waveformvalues_wholedt[iteration, 0, r] = getSource(time_x-start, freq, wavetype, dt)
|
||||
waveformvalues_wholedt[iteration, 1, r] = getSource(time_y-start, freq, wavetype, dt)
|
||||
waveformvalues_wholedt[iteration, 2, r] = getSource(time_z-start, freq, wavetype, dt)
|
||||
|
||||
|
||||
cpdef void updatePlaneWave(
|
||||
int n,
|
||||
int nthreads,
|
||||
float_or_double[:, ::1] H_fields,
|
||||
float_or_double[:, ::1] E_fields,
|
||||
float_or_double[:] updatecoeffsE,
|
||||
float_or_double[:] updatecoeffsH,
|
||||
float_or_double[:, :, ::1] Ex,
|
||||
float_or_double[:, :, ::1] Ey,
|
||||
float_or_double[:, :, ::1] Ez,
|
||||
float_or_double[:, :, ::1] Hx,
|
||||
float_or_double[:, :, ::1] Hy,
|
||||
float_or_double[:, :, ::1] Hz,
|
||||
double[:] projections,
|
||||
float_or_double[:, ::1] waveformvalues_wholedt,
|
||||
int[:] m,
|
||||
int[:] corners,
|
||||
bint precompute,
|
||||
int iteration,
|
||||
double dt,
|
||||
double ds,
|
||||
double c,
|
||||
double start,
|
||||
double stop,
|
||||
double freq,
|
||||
char* wavetype
|
||||
):
|
||||
initializeMagneticFields(m, H_fields, projections, waveformvalues_wholedt, precompute, iteration, dt, ds, c, start, stop, freq, wavetype)
|
||||
updateMagneticFields(n, H_fields, E_fields, updatecoeffsH, m)
|
||||
applyTFSFMagnetic(nthreads, Hx, Hy, Hz, E_fields, updatecoeffsH, m, corners)
|
||||
applyTFSFElectric(nthreads, Ex, Ey, Ez, H_fields, updatecoeffsE, m, corners)
|
||||
updateElectricFields(n, H_fields, E_fields, updatecoeffsE, m)
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cdef void takeSnapshot3D(double[:, :, ::1] field, char* filename):
|
||||
"""Writes fields of plane wave simulation at a particular time step.
|
||||
|
||||
Args:
|
||||
fields: double array of fields for grid cells over TFSF box at
|
||||
particular indices of TFSF box at particular time step.
|
||||
filename: string of file location where fields are to be written.
|
||||
|
||||
"""
|
||||
|
||||
cdef FILE *fptr = fopen(filename, "wb")
|
||||
fwrite(&field[0, 0, 0], sizeof(double), field.size, fptr)
|
||||
fclose(fptr)
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -47,7 +47,6 @@ def store_outputs(G: FDTDGrid, iteration: int):
|
||||
# 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:
|
||||
@@ -121,13 +120,21 @@ def write_hd5_data(basegrp, grid, is_subgrid=False):
|
||||
for srcindex, src in enumerate(srclist):
|
||||
grp = basegrp.create_group(f"srcs/src{str(srcindex + 1)}")
|
||||
grp.attrs["Type"] = type(src).__name__
|
||||
grp.attrs["Position"] = (src.xcoord * grid.dx, src.ycoord * grid.dy, src.zcoord * grid.dz)
|
||||
grp.attrs["Position"] = (
|
||||
src.xcoord * grid.dx,
|
||||
src.ycoord * grid.dy,
|
||||
src.zcoord * grid.dz,
|
||||
)
|
||||
|
||||
# Create group for transmission lines; add positional data, line resistance and
|
||||
# line discretisation attributes; write arrays for line voltages and currents
|
||||
for tlindex, tl in enumerate(grid.transmissionlines):
|
||||
grp = basegrp.create_group("tls/tl" + str(tlindex + 1))
|
||||
grp.attrs["Position"] = (tl.xcoord * grid.dx, tl.ycoord * grid.dy, tl.zcoord * grid.dz)
|
||||
grp.attrs["Position"] = (
|
||||
tl.xcoord * grid.dx,
|
||||
tl.ycoord * grid.dy,
|
||||
tl.zcoord * grid.dz,
|
||||
)
|
||||
grp.attrs["Resistance"] = tl.resistance
|
||||
grp.attrs["dl"] = tl.dl
|
||||
# Save incident voltage and current
|
||||
@@ -146,7 +153,62 @@ def write_hd5_data(basegrp, grid, is_subgrid=False):
|
||||
grp = basegrp.create_group("rxs/rx" + str(rxindex + 1))
|
||||
if rx.ID:
|
||||
grp.attrs["Name"] = rx.ID
|
||||
grp.attrs["Position"] = (rx.xcoord * grid.dx, rx.ycoord * grid.dy, rx.zcoord * grid.dz)
|
||||
grp.attrs["Position"] = (
|
||||
rx.xcoord * grid.dx,
|
||||
rx.ycoord * grid.dy,
|
||||
rx.zcoord * grid.dz,
|
||||
)
|
||||
|
||||
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
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -85,7 +85,10 @@ class FractalSurface:
|
||||
|
||||
# Positional vector at centre of array, scaled by weighting
|
||||
v1 = np.array(
|
||||
[self.weighting[0] * (surfacedims[0]) / 2, self.weighting[1] * (surfacedims[1]) / 2]
|
||||
[
|
||||
self.weighting[0] * (surfacedims[0]) / 2,
|
||||
self.weighting[1] * (surfacedims[1]) / 2,
|
||||
]
|
||||
)
|
||||
|
||||
# 2D array of random numbers to be convolved with the fractal function
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -36,7 +36,14 @@ from gprMax.fractals 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
|
||||
from gprMax.sources import HertzianDipole, MagneticDipole, Source, TransmissionLine, VoltageSource
|
||||
from gprMax.sources import (
|
||||
DiscretePlaneWave,
|
||||
HertzianDipole,
|
||||
MagneticDipole,
|
||||
Source,
|
||||
TransmissionLine,
|
||||
VoltageSource,
|
||||
)
|
||||
from gprMax.utilities.utilities import fft_power, get_terminal_width, round_value
|
||||
from gprMax.waveforms import Waveform
|
||||
|
||||
@@ -109,6 +116,7 @@ class FDTDGrid:
|
||||
self.hertziandipoles: List[HertzianDipole] = []
|
||||
self.magneticdipoles: List[MagneticDipole] = []
|
||||
self.transmissionlines: List[TransmissionLine] = []
|
||||
self.discreteplanewaves: List[DiscretePlaneWave] = []
|
||||
self.rxs: List[Rx] = []
|
||||
self.snapshots = [] # List[Snapshot]
|
||||
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -248,6 +248,7 @@ def check_cmd_names(processedlines, checkessential=True):
|
||||
"#hertzian_dipole",
|
||||
"#magnetic_dipole",
|
||||
"#transmission_line",
|
||||
"#discrete_plane_wave",
|
||||
"#excitation_file",
|
||||
"#rx",
|
||||
"#rx_array",
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -158,8 +158,7 @@ def process_geometrycmds(geometry):
|
||||
|
||||
# Isotropic case with user specified averaging
|
||||
elif len(tmp) == 9:
|
||||
averaging = check_averaging(tmp[8].lower())
|
||||
box = Box(p1=p1, p2=p2, material_id=tmp[7], averaging=averaging)
|
||||
box = Box(p1=p1, p2=p2, material_id=tmp[7], averaging=tmp[8].lower())
|
||||
|
||||
# Uniaxial anisotropic case
|
||||
elif len(tmp) == 10:
|
||||
@@ -186,8 +185,7 @@ def process_geometrycmds(geometry):
|
||||
|
||||
# Isotropic case with user specified averaging
|
||||
elif len(tmp) == 10:
|
||||
averaging = check_averaging(tmp[9].lower())
|
||||
cylinder = Cylinder(p1=p1, p2=p2, r=r, material_id=tmp[8], averaging=averaging)
|
||||
cylinder = Cylinder(p1=p1, p2=p2, r=r, material_id=tmp[8], averaging=tmp[9].lower())
|
||||
|
||||
# Uniaxial anisotropic case
|
||||
elif len(tmp) == 11:
|
||||
@@ -215,8 +213,14 @@ def process_geometrycmds(geometry):
|
||||
|
||||
# Isotropic case with user specified averaging
|
||||
elif len(tmp) == 11:
|
||||
averaging = check_averaging(tmp[10].lower())
|
||||
cone = Cone(p1=p1, p2=p2, r1=r1, r2=r2, material_id=tmp[9], averaging=averaging)
|
||||
cone = Cone(
|
||||
p1=p1,
|
||||
p2=p2,
|
||||
r1=r1,
|
||||
r2=r2,
|
||||
material_id=tmp[9],
|
||||
averaging=tmp[10].lower(),
|
||||
)
|
||||
|
||||
# Uniaxial anisotropic case
|
||||
elif len(tmp) == 12:
|
||||
@@ -258,7 +262,6 @@ def process_geometrycmds(geometry):
|
||||
|
||||
# Isotropic case with user specified averaging
|
||||
elif len(tmp) == 11:
|
||||
averaging = check_averaging(tmp[10].lower())
|
||||
cylindrical_sector = CylindricalSector(
|
||||
normal=normal,
|
||||
ctr1=ctr1,
|
||||
@@ -268,7 +271,7 @@ def process_geometrycmds(geometry):
|
||||
r=r,
|
||||
start=start,
|
||||
end=end,
|
||||
averaging=averaging,
|
||||
averaging=tmp[10].lower(),
|
||||
material_id=tmp[9],
|
||||
)
|
||||
|
||||
@@ -306,8 +309,7 @@ def process_geometrycmds(geometry):
|
||||
|
||||
# Isotropic case with user specified averaging
|
||||
elif len(tmp) == 7:
|
||||
averaging = check_averaging(tmp[6].lower())
|
||||
sphere = Sphere(p1=p1, r=r, material_id=tmp[5], averaging=averaging)
|
||||
sphere = Sphere(p1=p1, r=r, material_id=tmp[5], averaging=tmp[6].lower())
|
||||
|
||||
# Uniaxial anisotropic case
|
||||
elif len(tmp) == 8:
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -23,6 +23,7 @@ from .user_objects.cmds_multiuse import (
|
||||
AddDebyeDispersion,
|
||||
AddDrudeDispersion,
|
||||
AddLorentzDispersion,
|
||||
DiscretePlaneWave,
|
||||
ExcitationFile,
|
||||
HertzianDipole,
|
||||
MagneticDipole,
|
||||
@@ -201,6 +202,55 @@ def process_multicmds(multicmds):
|
||||
|
||||
scene_objects.append(tl)
|
||||
|
||||
cmdname = "#discrete_plane_wave"
|
||||
if multicmds[cmdname] is not None:
|
||||
for cmdinstance in multicmds[cmdname]:
|
||||
tmp = cmdinstance.split()
|
||||
if len(tmp) < 10:
|
||||
logger.exception(
|
||||
"'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least ten parameters"
|
||||
)
|
||||
raise ValueError
|
||||
|
||||
if len(tmp) == 10:
|
||||
plWave = DiscretePlaneWave(
|
||||
p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])),
|
||||
p2=(float(tmp[3]), float(tmp[4]), float(tmp[5])),
|
||||
psi=float(tmp[6]),
|
||||
phi=float(tmp[7]),
|
||||
theta=float(tmp[8]),
|
||||
waveform_id=tmp[9],
|
||||
)
|
||||
elif len(tmp) == 11:
|
||||
plWave = DiscretePlaneWave(
|
||||
p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])),
|
||||
p2=(float(tmp[3]), float(tmp[4]), float(tmp[5])),
|
||||
psi=float(tmp[6]),
|
||||
phi=float(tmp[7]),
|
||||
theta=float(tmp[8]),
|
||||
waveform_id=tmp[9],
|
||||
material_ID=int(tmp[10]),
|
||||
)
|
||||
elif len(tmp) == 13:
|
||||
plWave = DiscretePlaneWave(
|
||||
p1=(float(tmp[0]), float(tmp[1]), float(tmp[2])),
|
||||
p2=(float(tmp[3]), float(tmp[4]), float(tmp[5])),
|
||||
psi=float(tmp[6]),
|
||||
phi=float(tmp[7]),
|
||||
theta=float(tmp[8]),
|
||||
waveform_id=tmp[9],
|
||||
material_ID=int(tmp[10]),
|
||||
start=float(tmp[11]),
|
||||
stop=float(tmp[12]),
|
||||
)
|
||||
else:
|
||||
logger.exception(
|
||||
"'" + cmdname + ": " + " ".join(tmp) + "'" + " too many parameters"
|
||||
)
|
||||
raise ValueError
|
||||
|
||||
scene_objects.append(plWave)
|
||||
|
||||
cmdname = "#excitation_file"
|
||||
if multicmds[cmdname] is not None:
|
||||
for cmdinstance in multicmds[cmdname]:
|
||||
@@ -366,7 +416,11 @@ def process_multicmds(multicmds):
|
||||
alpha.append(float(tmp[pole + 2]))
|
||||
|
||||
lorentz_dispersion = AddLorentzDispersion(
|
||||
poles=poles, material_ids=material_ids, er_delta=er_delta, tau=tau, alpha=alpha
|
||||
poles=poles,
|
||||
material_ids=material_ids,
|
||||
er_delta=er_delta,
|
||||
tau=tau,
|
||||
alpha=alpha,
|
||||
)
|
||||
scene_objects.append(lorentz_dispersion)
|
||||
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -252,7 +252,7 @@ class PML:
|
||||
kappamin = sum(cfs.kappa.min for cfs in self.CFS)
|
||||
if kappamin < 1:
|
||||
logger.exception(
|
||||
f"Sum of kappamin value(s) for PML is {kappamin} " "and must be greater than one."
|
||||
f"Sum of kappamin value(s) for PML is {kappamin} and must be greater than one."
|
||||
)
|
||||
raise ValueError
|
||||
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -65,6 +65,7 @@ class Solver:
|
||||
self.updates.update_magnetic()
|
||||
self.updates.update_magnetic_pml()
|
||||
self.updates.update_magnetic_sources(iteration)
|
||||
# self.updates.update_plane_waves()
|
||||
if isinstance(self.updates, MPIUpdates):
|
||||
self.updates.halo_swap_magnetic()
|
||||
if isinstance(self.updates, SubgridUpdates):
|
||||
|
||||
698
gprMax/sources.py
普通文件 -> 可执行文件
698
gprMax/sources.py
普通文件 -> 可执行文件
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -24,6 +24,13 @@ import numpy.typing as npt
|
||||
import gprMax.config as config
|
||||
from gprMax.waveforms import Waveform
|
||||
|
||||
from .cython.plane_wave import (
|
||||
calculate1DWaveformValues,
|
||||
getIntegerForAngles,
|
||||
getProjections,
|
||||
getSource,
|
||||
updatePlaneWave,
|
||||
)
|
||||
from .utilities.utilities import round_value
|
||||
|
||||
|
||||
@@ -37,7 +44,11 @@ class Source:
|
||||
self.coordorigin = np.zeros(3, dtype=np.int32)
|
||||
self.start = 0.0
|
||||
self.stop = 0.0
|
||||
self.waveform: Waveform
|
||||
self.waveformID = None
|
||||
# Waveform values for sources that need to be calculated on whole timesteps
|
||||
self.waveformvalues_wholedt = None
|
||||
# Waveform values for sources that need to be calculated on half timesteps
|
||||
self.waveformvalues_halfdt = None
|
||||
|
||||
@property
|
||||
def xcoord(self) -> int:
|
||||
@@ -87,33 +98,6 @@ class Source:
|
||||
def zcoordorigin(self, value: int):
|
||||
self.coordorigin[2] = value
|
||||
|
||||
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(
|
||||
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(
|
||||
iterations, dtype=config.sim_config.dtypes["float_or_double"]
|
||||
)
|
||||
|
||||
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] = self.waveform.calculate_value(time, dt)
|
||||
self.waveformvalues_halfdt[iteration] = self.waveform.calculate_value(
|
||||
time + 0.5 * dt, dt
|
||||
)
|
||||
|
||||
|
||||
class VoltageSource(Source):
|
||||
"""A voltage source can be a hard source if it's resistance is zero,
|
||||
@@ -126,6 +110,45 @@ class VoltageSource(Source):
|
||||
super().__init__()
|
||||
self.resistance = None
|
||||
|
||||
def calculate_waveform_values(self, G):
|
||||
"""Calculates all waveform values for source for duration of simulation.
|
||||
|
||||
Args:
|
||||
G: FDTDGrid class describing a grid in a model.
|
||||
"""
|
||||
|
||||
# Check if a source matches existing source in terms of waveform and
|
||||
# does not have a customised start/stop time. If so, use its
|
||||
# pre-calculated waveform values, otherwise calculate them.
|
||||
src_match = False
|
||||
|
||||
if self.start == 0 and self.stop == G.timewindow:
|
||||
for src in G.voltagesources:
|
||||
if src.waveformID == self.waveformID:
|
||||
src_match = True
|
||||
self.waveformvalues_halfdt = src.waveformvalues_halfdt
|
||||
self.waveformvalues_wholedt = src.waveformvalues_wholedt
|
||||
|
||||
if not src_match:
|
||||
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
|
||||
self.waveformvalues_halfdt = np.zeros(
|
||||
(G.iterations), dtype=config.sim_config.dtypes["float_or_double"]
|
||||
)
|
||||
self.waveformvalues_wholedt = np.zeros(
|
||||
(G.iterations), dtype=config.sim_config.dtypes["float_or_double"]
|
||||
)
|
||||
|
||||
for iteration in range(G.iterations):
|
||||
time = G.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_halfdt[iteration] = waveform.calculate_value(
|
||||
time + 0.5 * G.dt, G.dt
|
||||
)
|
||||
self.waveformvalues_wholedt[iteration] = waveform.calculate_value(time, G.dt)
|
||||
|
||||
def update_electric(self, iteration, updatecoeffsE, ID, Ex, Ey, Ez, G):
|
||||
"""Updates electric field values for a voltage source.
|
||||
|
||||
@@ -217,6 +240,40 @@ class HertzianDipole(Source):
|
||||
super().__init__()
|
||||
self.dl = 0.0
|
||||
|
||||
def calculate_waveform_values(self, G):
|
||||
"""Calculates all waveform values for source for duration of simulation.
|
||||
|
||||
Args:
|
||||
G: FDTDGrid class describing a grid in a model.
|
||||
"""
|
||||
|
||||
# Check if a source matches existing source in terms of waveform and
|
||||
# does not have a customised start/stop time. If so, use its
|
||||
# pre-calculated waveform values, otherwise calculate them.
|
||||
src_match = False
|
||||
|
||||
if self.start == 0 and self.stop == G.timewindow:
|
||||
for src in G.hertziandipoles:
|
||||
if src.waveformID == self.waveformID:
|
||||
src_match = True
|
||||
self.waveformvalues_halfdt = src.waveformvalues_halfdt
|
||||
|
||||
if not src_match:
|
||||
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
|
||||
self.waveformvalues_halfdt = np.zeros(
|
||||
(G.iterations), dtype=config.sim_config.dtypes["float_or_double"]
|
||||
)
|
||||
|
||||
for iteration in range(G.iterations):
|
||||
time = G.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_halfdt[iteration] = waveform.calculate_value(
|
||||
time + 0.5 * G.dt, G.dt
|
||||
)
|
||||
|
||||
def update_electric(self, iteration, updatecoeffsE, ID, Ex, Ey, Ez, G):
|
||||
"""Updates electric field values for a Hertzian dipole.
|
||||
|
||||
@@ -263,6 +320,38 @@ class HertzianDipole(Source):
|
||||
class MagneticDipole(Source):
|
||||
"""A magnetic dipole is an additive source (magnetic current density)."""
|
||||
|
||||
def calculate_waveform_values(self, G):
|
||||
"""Calculates all waveform values for source for duration of simulation.
|
||||
|
||||
Args:
|
||||
G: FDTDGrid class describing a grid in a model.
|
||||
"""
|
||||
|
||||
# Check if a source matches existing source in terms of waveform and
|
||||
# does not have a customised start/stop time. If so, use its
|
||||
# pre-calculated waveform values, otherwise calculate them.
|
||||
src_match = False
|
||||
|
||||
if self.start == 0 and self.stop == G.timewindow:
|
||||
for src in G.magneticdipoles:
|
||||
if src.waveformID == self.waveformID:
|
||||
src_match = True
|
||||
self.waveformvalues_wholedt = src.waveformvalues_wholedt
|
||||
|
||||
if not src_match:
|
||||
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
|
||||
self.waveformvalues_wholedt = np.zeros(
|
||||
(G.iterations), dtype=config.sim_config.dtypes["float_or_double"]
|
||||
)
|
||||
|
||||
for iteration in range(G.iterations):
|
||||
time = G.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)
|
||||
|
||||
def update_magnetic(self, iteration, updatecoeffsH, ID, Hx, Hy, Hz, G):
|
||||
"""Updates magnetic field values for a magnetic dipole.
|
||||
|
||||
@@ -411,6 +500,45 @@ class TransmissionLine(Source):
|
||||
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_waveform_values(self, G):
|
||||
"""Calculates all waveform values for source for duration of simulation.
|
||||
|
||||
Args:
|
||||
G: FDTDGrid class describing a grid in a model.
|
||||
"""
|
||||
|
||||
# Check if a source matches existing source in terms of waveform and
|
||||
# does not have a customised start/stop time. If so, use its
|
||||
# pre-calculated waveform values, otherwise calculate them.
|
||||
src_match = False
|
||||
|
||||
if self.start == 0 and self.stop == G.timewindow:
|
||||
for src in G.transmissionlines:
|
||||
if src.waveformID == self.waveformID:
|
||||
src_match = True
|
||||
self.waveformvalues_wholedt = src.waveformvalues_wholedt
|
||||
self.waveformvalues_halfdt = src.waveformvalues_halfdt
|
||||
|
||||
if not src_match:
|
||||
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
|
||||
self.waveformvalues_wholedt = np.zeros(
|
||||
(G.iterations), dtype=config.sim_config.dtypes["float_or_double"]
|
||||
)
|
||||
self.waveformvalues_halfdt = np.zeros(
|
||||
(G.iterations), dtype=config.sim_config.dtypes["float_or_double"]
|
||||
)
|
||||
|
||||
for iteration in range(G.iterations):
|
||||
time = G.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
|
||||
)
|
||||
|
||||
def calculate_incident_V_I(self, G):
|
||||
"""Calculates the incident voltage and current with a long length
|
||||
transmission line not connected to the main grid
|
||||
@@ -547,3 +675,515 @@ class TransmissionLine(Source):
|
||||
self.current[self.antpos] = G.calculate_Iz(i, j, k)
|
||||
|
||||
self.update_current(iteration, G)
|
||||
|
||||
|
||||
class DiscretePlaneWave(Source):
|
||||
"""Implements the discrete plane wave (DPW) formulation as described in
|
||||
Tan, T.; Potter, M. (2010). FDTD Discrete Planewave (FDTD-DPW)
|
||||
Formulation for a Perfectly Matched Source in TFSF Simulations., 58(8),
|
||||
0–2648. doi:10.1109/tap.2010.2050446
|
||||
"""
|
||||
|
||||
def __init__(self, G):
|
||||
"""
|
||||
Args:
|
||||
m: int array stores the integer mappings, m_x, m_y, m_z which
|
||||
determine the rational angles last element stores
|
||||
max(m_x, m_y, m_z).
|
||||
directions: int array stores the directions of propagation of DPW.
|
||||
dimensions: int stores the number of dimensions in which the
|
||||
simulation is run (2D or 3D).
|
||||
time_dimension: int stores the time length over which the simulation
|
||||
is run.
|
||||
E_fields: double array stores the electric flieds associated with
|
||||
1D DPW.
|
||||
H_fields: double array stores the magnetic fields associated with
|
||||
1D DPW.
|
||||
G: FDTDGrid class describing a grid in a model.
|
||||
"""
|
||||
|
||||
super().__init__()
|
||||
self.m = np.zeros(3 + 1, dtype=np.int32) # +1 to store the max(m_x, m_y, m_z)
|
||||
self.directions = np.zeros(3, dtype=np.int32)
|
||||
self.length = 0
|
||||
self.projections = np.zeros(3, dtype=config.sim_config.dtypes["float_or_double"])
|
||||
self.corners = None
|
||||
self.materialID = 1
|
||||
self.ds = 0
|
||||
|
||||
def initializeDiscretePlaneWave(self, psi, phi, Delta_phi, theta, Delta_theta, G):
|
||||
"""Creates a DPW, assigns memory to the grids, and gets field values
|
||||
at different time and space indices.
|
||||
|
||||
Args:
|
||||
psi: float for polarization angle of the incident plane wave.
|
||||
phi: float for azimuthal angle (radians) of the incident plane wave.
|
||||
Delta_phi: float for permissible error in the rational angle
|
||||
(radians) approximation to phi.
|
||||
theta: float for polar angle (radians) of the incident plane wave.
|
||||
Delta_theta: float for permissible error in the rational angle
|
||||
(radians) approximation to theta.
|
||||
G: FDTDGrid class describing a grid in a model.
|
||||
number: int for number of cells in the 3D FDTD simulation.
|
||||
dx: double for separation between adjacent cells in the x direction.
|
||||
dy: double for separation between adjacent cells in the y direction.
|
||||
dz: double for separation between adjacent cells in the z direction.
|
||||
dt: double for time step for the FDTD simulation.
|
||||
|
||||
Returns:
|
||||
E_fields: double array for electric field for the DPW as it evolves
|
||||
over time and space indices.
|
||||
H_fields: double array for magnetic field for the DPW as it evolves
|
||||
over time and space indices.
|
||||
C: double array stores coefficients of the fields for the update
|
||||
equation of the electric fields.
|
||||
D: double array stores coefficients of the fields for the update
|
||||
equation of the magnetic fields.
|
||||
"""
|
||||
|
||||
# Get the integers for the nearest rational angle
|
||||
self.directions, self.m[:3] = getIntegerForAngles(
|
||||
phi, Delta_phi, theta, Delta_theta, np.array([G.dx, G.dy, G.dz])
|
||||
)
|
||||
|
||||
# Store max(m_x, m_y, m_z) in the last element of the array
|
||||
self.m[3] = max(self.m)
|
||||
|
||||
# Set an appropriate length fo the one dimensional arrays
|
||||
self.length = int(2 * max(self.m[:-1]) * G.iterations)
|
||||
|
||||
self.E_fields = np.zeros(
|
||||
(3, self.length),
|
||||
order="C",
|
||||
dtype=config.sim_config.dtypes["float_or_double"],
|
||||
)
|
||||
self.H_fields = np.zeros(
|
||||
(3, self.length),
|
||||
order="C",
|
||||
dtype=config.sim_config.dtypes["float_or_double"],
|
||||
)
|
||||
|
||||
# 1D grid has no ABC to terminate it, sufficiently long array prevents
|
||||
# reflections from the back
|
||||
# Projections for field components
|
||||
projections_h, P = getProjections(psi * 180 / np.pi, self.m)
|
||||
# Scale the projection vector for the magnetic field
|
||||
self.projections = projections_h / np.sqrt(config.m0 / config.e0)
|
||||
|
||||
if self.m[0] == 0: # Calculate dr that is needed for sourcing the 1D array
|
||||
if self.m[1] == 0:
|
||||
if self.m[2] == 0:
|
||||
raise ValueError("not all m_i values can be zero")
|
||||
else:
|
||||
self.ds = P[2] * G.dz / self.m[2]
|
||||
else:
|
||||
self.ds = P[1] * G.dy / self.m[1]
|
||||
else:
|
||||
self.ds = P[0] * G.dx / self.m[0]
|
||||
|
||||
def calculate_waveform_values(self, G, cythonize=True):
|
||||
"""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, 3, self.m[3]),
|
||||
dtype=config.sim_config.dtypes["float_or_double"],
|
||||
)
|
||||
|
||||
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
|
||||
if cythonize:
|
||||
calculate1DWaveformValues(
|
||||
self.waveformvalues_wholedt,
|
||||
G.iterations,
|
||||
self.m,
|
||||
G.dt,
|
||||
self.ds,
|
||||
config.c,
|
||||
self.start,
|
||||
self.stop,
|
||||
waveform.freq,
|
||||
waveform.type.encode("UTF-8"),
|
||||
)
|
||||
else:
|
||||
for dimension in range(3):
|
||||
for iteration in range(G.iterations):
|
||||
for r in range(self.m[3]):
|
||||
time = (
|
||||
G.dt * iteration
|
||||
- (
|
||||
r
|
||||
+ (self.m[(dimension + 1) % 3] + self.m[(dimension + 2) % 3]) * 0.5
|
||||
)
|
||||
* self.ds
|
||||
/ config.c
|
||||
)
|
||||
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, dimension, r
|
||||
] = waveform.calculate_value(time, G.dt)
|
||||
|
||||
def update_plane_wave(
|
||||
self,
|
||||
nthreads,
|
||||
updatecoeffsE,
|
||||
updatecoeffsH,
|
||||
Ex,
|
||||
Ey,
|
||||
Ez,
|
||||
Hx,
|
||||
Hy,
|
||||
Hz,
|
||||
iteration,
|
||||
G,
|
||||
cythonize=True,
|
||||
precompute=True,
|
||||
):
|
||||
if cythonize:
|
||||
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
|
||||
updatePlaneWave(
|
||||
self.length,
|
||||
nthreads,
|
||||
self.H_fields,
|
||||
self.E_fields,
|
||||
updatecoeffsE[self.material_ID, :],
|
||||
updatecoeffsH[self.material_ID, :],
|
||||
Ex,
|
||||
Ey,
|
||||
Ez,
|
||||
Hx,
|
||||
Hy,
|
||||
Hz,
|
||||
self.projections,
|
||||
self.waveformvalues_wholedt[iteration, :, :],
|
||||
self.m,
|
||||
self.corners,
|
||||
precompute,
|
||||
iteration,
|
||||
G.dt,
|
||||
self.ds,
|
||||
config.c,
|
||||
self.start,
|
||||
self.stop,
|
||||
waveform.freq,
|
||||
waveform.type.encode("UTF-8"),
|
||||
)
|
||||
else:
|
||||
self.update_magnetic_field_1D(G, precompute)
|
||||
self.apply_TFSF_conditions_magnetic(G)
|
||||
self.apply_TFSF_conditions_electric(G)
|
||||
self.update_electric_field_1D(G)
|
||||
|
||||
def initialize_magnetic_fields_1D(self, G, precompute):
|
||||
if precompute:
|
||||
for dimension in range(3):
|
||||
for r in range(self.m[3]):
|
||||
# Assign source values of magnetic field to first few gridpoints
|
||||
self.H_fields[dimension, r] = (
|
||||
self.projections[dimension]
|
||||
* self.waveformvalues_wholedt[G.iteration, dimension, r]
|
||||
)
|
||||
# self.getSource(self.real_time - (j+(self.m[(i+1)%3]+self.m[(i+2)%3])*0.5)*self.ds/config.c)#, self.waveformID, G.dt)
|
||||
else:
|
||||
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
|
||||
for dimension in range(3):
|
||||
for r in range(self.m[3]):
|
||||
# Assign source values of magnetic field to first few gridpoints
|
||||
self.H_fields[dimension, r] = self.projections[dimension] * getSource(
|
||||
G.iteration * G.dt
|
||||
- (r + (self.m[(dimension + 1) % 3] + self.m[(dimension + 2) % 3]) * 0.5)
|
||||
* self.ds
|
||||
/ config.c,
|
||||
waveform.freq,
|
||||
waveform.type.encode("UTF-8"),
|
||||
G.dt,
|
||||
)
|
||||
|
||||
def update_magnetic_field_1D(self, G, precompute=True):
|
||||
"""Updates magnetic fields for the next time step using Equation 8 of
|
||||
DOI: 10.1109/LAWP.2009.2016851
|
||||
|
||||
Args:
|
||||
n: int stores spatial length of the DPW array so that each length
|
||||
grid cell is updated when updateMagneticFields() called.
|
||||
H_coefficients: double array stores coefficients of the fields in
|
||||
the update equation for the magnetic field.
|
||||
H_fields: double array stores magnetic fields of the DPW until
|
||||
temporal index time.
|
||||
E_fields: double array stores electric fields of the DPW until
|
||||
temporal index time.
|
||||
time: int time index storing current axis number which would be
|
||||
updated for the H_fields.
|
||||
|
||||
Returns:
|
||||
H_fields: double array for magnetic field with the axis entry for
|
||||
the current time added.
|
||||
"""
|
||||
|
||||
self.initialize_magnetic_fields_1D(G, precompute)
|
||||
|
||||
for i in range(3): # Update each component of magnetic field
|
||||
materialH = G.ID[
|
||||
3 + i,
|
||||
(self.corners[0] + self.corners[3]) // 2,
|
||||
(self.corners[1] + self.corners[4]) // 2,
|
||||
(self.corners[2] + self.corners[5]) // 2,
|
||||
]
|
||||
# Update magnetic field at each spatial index
|
||||
for j in range(self.m[-1], self.length - self.m[-1]):
|
||||
self.H_fields[i, j] = (
|
||||
G.updatecoeffsH[materialH, 0] * self.H_fields[i, j]
|
||||
+ G.updatecoeffsH[materialH, (i + 2) % 3 + 1]
|
||||
* (
|
||||
self.E_fields[(i + 1) % 3, j + self.m[(i + 2) % 3]]
|
||||
- self.E_fields[(i + 1) % 3, j]
|
||||
)
|
||||
- G.updatecoeffsH[materialH, (i + 1) % 3 + 1]
|
||||
* (
|
||||
self.E_fields[(i + 2) % 3, j + self.m[(i + 1) % 3]]
|
||||
- self.E_fields[(i + 2) % 3, j]
|
||||
)
|
||||
) # equation 8 of Tan, Potter paper
|
||||
|
||||
def update_electric_field_1D(self, G):
|
||||
"""Updates electric fields for the next time step using Equation 9 of
|
||||
DOI: 10.1109/LAWP.2009.2016851
|
||||
|
||||
Args:
|
||||
n: int stores spatial length of DPW array so that each length grid
|
||||
cell is updated when updateMagneticFields() is called.
|
||||
E_coefficients: double array stores coefficients of the fields in
|
||||
the update equation for the electric field.
|
||||
H_fields: double array stores magnetic fields of the DPW until
|
||||
temporal index time.
|
||||
E_fields: double array stores electric fields of the DPW until
|
||||
temporal index time.
|
||||
time: int time index storing current axis number which would be
|
||||
updated for the E_fields.
|
||||
|
||||
Returns:
|
||||
E_fields: double array for electric field with the axis entry for
|
||||
the current time added.
|
||||
|
||||
"""
|
||||
|
||||
for i in range(3): # Update each component of electric field
|
||||
materialE = G.ID[
|
||||
i,
|
||||
(self.corners[0] + self.corners[3]) // 2,
|
||||
(self.corners[1] + self.corners[4]) // 2,
|
||||
(self.corners[2] + self.corners[5]) // 2,
|
||||
]
|
||||
# Update electric field at each spatial index
|
||||
for j in range(self.m[-1], self.length - self.m[-1]):
|
||||
self.E_fields[i, j] = (
|
||||
G.updatecoeffsE[materialE, 0] * self.E_fields[i, j]
|
||||
+ G.updatecoeffsE[materialE, (i + 2) % 3 + 1]
|
||||
* (
|
||||
self.H_fields[(i + 2) % 3, j]
|
||||
- self.H_fields[(i + 2) % 3, j - self.m[(i + 1) % 3]]
|
||||
)
|
||||
- G.updatecoeffsE[materialE, (i + 1) % 3 + 1]
|
||||
* (
|
||||
self.H_fields[(i + 1) % 3, j]
|
||||
- self.H_fields[(i + 1) % 3, j - self.m[(i + 2) % 3]]
|
||||
)
|
||||
) # equation 9 of Tan, Potter paper
|
||||
|
||||
def getField(self, i, j, k, array, m, component):
|
||||
return array[component, np.dot(m[:-1], np.array([i, j, k]))]
|
||||
|
||||
def apply_TFSF_conditions_magnetic(self, G):
|
||||
# **** constant x faces -- scattered-field nodes ****
|
||||
i = self.corners[0]
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
# correct Hy at firstX-1/2 by subtracting Ez_inc
|
||||
G.Hy[i - 1, j, k] -= G.updatecoeffsH[G.ID[4, i, j, k], 1] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 2
|
||||
)
|
||||
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Hz at firstX-1/2 by adding Ey_inc
|
||||
G.Hz[i - 1, j, k] += G.updatecoeffsH[G.ID[5, i, j, k], 1] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 1
|
||||
)
|
||||
|
||||
i = self.corners[3]
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
# correct Hy at lastX+1/2 by adding Ez_inc
|
||||
G.Hy[i, j, k] += G.updatecoeffsH[G.ID[4, i, j, k], 1] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 2
|
||||
)
|
||||
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Hz at lastX+1/2 by subtractinging Ey_inc
|
||||
G.Hz[i, j, k] -= G.updatecoeffsH[G.ID[5, i, j, k], 1] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 1
|
||||
)
|
||||
|
||||
# **** constant y faces -- scattered-field nodes ****
|
||||
j = self.corners[1]
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
# correct Hx at firstY-1/2 by adding Ez_inc
|
||||
G.Hx[i, j - 1, k] += G.updatecoeffsH[G.ID[3, i, j, k], 2] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 2
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Hz at firstY-1/2 by subtracting Ex_inc
|
||||
G.Hz[i, j - 1, k] -= G.updatecoeffsH[G.ID[5, i, j, k], 2] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 0
|
||||
)
|
||||
|
||||
j = self.corners[4]
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
# correct Hx at lastY+1/2 by subtracting Ez_inc
|
||||
G.Hx[i, j, k] -= G.updatecoeffsH[G.ID[3, i, j, k], 2] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 2
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Hz at lastY-1/2 by adding Ex_inc
|
||||
G.Hz[i, j, k] += G.updatecoeffsH[G.ID[5, i, j, k], 2] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 0
|
||||
)
|
||||
|
||||
# **** constant z faces -- scattered-field nodes ****
|
||||
k = self.corners[2]
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
# correct Hy at firstZ-1/2 by adding Ex_inc
|
||||
G.Hy[i, j, k - 1] += G.updatecoeffsH[G.ID[4, i, j, k], 3] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 0
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
# correct Hx at firstZ-1/2 by subtracting Ey_inc
|
||||
G.Hx[i, j, k - 1] -= G.updatecoeffsH[G.ID[3, i, j, k], 3] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 1
|
||||
)
|
||||
|
||||
k = self.corners[5]
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
# correct Hy at firstZ-1/2 by subtracting Ex_inc
|
||||
G.Hy[i, j, k] -= G.updatecoeffsH[G.ID[4, i, j, k], 3] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 0
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
# correct Hx at lastZ+1/2 by adding Ey_inc
|
||||
G.Hx[i, j, k] += G.updatecoeffsH[G.ID[3, i, j, k], 3] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 1
|
||||
)
|
||||
|
||||
def apply_TFSF_conditions_electric(self, G):
|
||||
# **** constant x faces -- total-field nodes ****/
|
||||
i = self.corners[0]
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
# correct Ez at firstX face by subtracting Hy_inc
|
||||
G.Ez[i, j, k] -= G.updatecoeffsE[G.ID[2, i, j, k], 1] * self.getField(
|
||||
i - 1, j, k, self.H_fields, self.m, 1
|
||||
)
|
||||
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Ey at firstX face by adding Hz_inc
|
||||
G.Ey[i, j, k] += G.updatecoeffsE[G.ID[1, i, j, k], 1] * self.getField(
|
||||
i - 1, j, k, self.H_fields, self.m, 2
|
||||
)
|
||||
|
||||
i = self.corners[3]
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
# correct Ez at lastX face by adding Hy_inc
|
||||
G.Ez[i, j, k] += G.updatecoeffsE[G.ID[2, i, j, k], 1] * self.getField(
|
||||
i, j, k, self.H_fields, self.m, 1
|
||||
)
|
||||
|
||||
i = self.corners[3]
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Ey at lastX face by subtracting Hz_inc
|
||||
G.Ey[i, j, k] -= G.updatecoeffsE[G.ID[1, i, j, k], 1] * self.getField(
|
||||
i, j, k, self.H_fields, self.m, 2
|
||||
)
|
||||
|
||||
# **** constant y faces -- total-field nodes ****/
|
||||
j = self.corners[1]
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
# correct Ez at firstY face by adding Hx_inc
|
||||
G.Ez[i, j, k] += G.updatecoeffsE[G.ID[2, i, j, k], 2] * self.getField(
|
||||
i, j - 1, k, self.H_fields, self.m, 0
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Ex at firstY face by subtracting Hz_inc
|
||||
G.Ex[i, j, k] -= G.updatecoeffsE[G.ID[0, i, j, k], 2] * self.getField(
|
||||
i, j - 1, k, self.H_fields, self.m, 2
|
||||
)
|
||||
|
||||
j = self.corners[4]
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
# correct Ez at lastY face by subtracting Hx_inc
|
||||
G.Ez[i, j, k] -= G.updatecoeffsE[G.ID[2, i, j, k], 2] * self.getField(
|
||||
i, j, k, self.H_fields, self.m, 0
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Ex at lastY face by adding Hz_inc
|
||||
G.Ex[i, j, k] += G.updatecoeffsE[G.ID[0, i, j, k], 2] * self.getField(
|
||||
i, j, k, self.H_fields, self.m, 2
|
||||
)
|
||||
|
||||
# **** constant z faces -- total-field nodes ****/
|
||||
k = self.corners[2]
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
# correct Ey at firstZ face by subtracting Hx_inc
|
||||
G.Ey[i, j, k] -= G.updatecoeffsE[G.ID[1, i, j, k], 3] * self.getField(
|
||||
i, j, k - 1, self.H_fields, self.m, 0
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
# correct Ex at firstZ face by adding Hy_inc
|
||||
G.Ex[i, j, k] += G.updatecoeffsE[G.ID[0, i, j, k], 3] * self.getField(
|
||||
i, j, k - 1, self.H_fields, self.m, 1
|
||||
)
|
||||
|
||||
k = self.corners[5]
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
# correct Ey at lastZ face by adding Hx_inc
|
||||
G.Ey[i, j, k] += G.updatecoeffsE[G.ID[1, i, j, k], 3] * self.getField(
|
||||
i, j, k, self.H_fields, self.m, 0
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
# correct Ex at lastZ face by subtracting Hy_inc
|
||||
G.Ex[i, j, k] -= G.updatecoeffsE[G.ID[0, i, j, k], 3] * self.getField(
|
||||
i, j, k, self.H_fields, self.m, 1
|
||||
)
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -233,7 +233,9 @@ class PrecursorNodesBase:
|
||||
|
||||
def interpolate_to_sub_grid(self, field, coords):
|
||||
x, z, x_sg, z_sg = coords
|
||||
interp_f = interpolate.RectBivariateSpline(x, z, field, kx=self.interpolation, ky=self.interpolation)
|
||||
interp_f = interpolate.RectBivariateSpline(
|
||||
x, z, field, kx=self.interpolation, ky=self.interpolation
|
||||
)
|
||||
f_i = interp_f(x_sg, z_sg)
|
||||
return f_i
|
||||
|
||||
@@ -417,18 +419,78 @@ class PrecursorNodes(PrecursorNodesBase):
|
||||
|
||||
# Spatially interpolate nodes
|
||||
slices = [
|
||||
["ex_front_1", True, (slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)), self.Ex],
|
||||
["ex_back_1", True, (slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)), self.Ex],
|
||||
["ez_front_1", False, (slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)), self.Ez],
|
||||
["ez_back_1", False, (slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)), self.Ez],
|
||||
["ey_left_1", True, (self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)), self.Ey],
|
||||
["ey_right_1", True, (self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)), self.Ey],
|
||||
["ez_left_1", False, (self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)), self.Ez],
|
||||
["ez_right_1", False, (self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)), self.Ez],
|
||||
["ex_bottom_1", True, (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0), self.Ex],
|
||||
["ex_top_1", True, (slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1), self.Ex],
|
||||
["ey_bottom_1", False, (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0), self.Ey],
|
||||
["ey_top_1", False, (slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1), self.Ey],
|
||||
[
|
||||
"ex_front_1",
|
||||
True,
|
||||
(slice(i0, i1, 1), self.j0, slice(k0, k1 + 1, 1)),
|
||||
self.Ex,
|
||||
],
|
||||
[
|
||||
"ex_back_1",
|
||||
True,
|
||||
(slice(i0, i1, 1), self.j1, slice(k0, k1 + 1, 1)),
|
||||
self.Ex,
|
||||
],
|
||||
[
|
||||
"ez_front_1",
|
||||
False,
|
||||
(slice(i0, i1 + 1, 1), self.j0, slice(k0, k1, 1)),
|
||||
self.Ez,
|
||||
],
|
||||
[
|
||||
"ez_back_1",
|
||||
False,
|
||||
(slice(i0, i1 + 1, 1), self.j1, slice(k0, k1, 1)),
|
||||
self.Ez,
|
||||
],
|
||||
[
|
||||
"ey_left_1",
|
||||
True,
|
||||
(self.i0, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
|
||||
self.Ey,
|
||||
],
|
||||
[
|
||||
"ey_right_1",
|
||||
True,
|
||||
(self.i1, slice(j0, j1, 1), slice(k0, k1 + 1, 1)),
|
||||
self.Ey,
|
||||
],
|
||||
[
|
||||
"ez_left_1",
|
||||
False,
|
||||
(self.i0, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
|
||||
self.Ez,
|
||||
],
|
||||
[
|
||||
"ez_right_1",
|
||||
False,
|
||||
(self.i1, slice(j0, j1 + 1, 1), slice(k0, k1, 1)),
|
||||
self.Ez,
|
||||
],
|
||||
[
|
||||
"ex_bottom_1",
|
||||
True,
|
||||
(slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k0),
|
||||
self.Ex,
|
||||
],
|
||||
[
|
||||
"ex_top_1",
|
||||
True,
|
||||
(slice(i0, i1, 1), slice(j0, j1 + 1, 1), self.k1),
|
||||
self.Ex,
|
||||
],
|
||||
[
|
||||
"ey_bottom_1",
|
||||
False,
|
||||
(slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k0),
|
||||
self.Ey,
|
||||
],
|
||||
[
|
||||
"ey_top_1",
|
||||
False,
|
||||
(slice(i0, i1 + 1, 1), slice(j0, j1, 1), self.k1),
|
||||
self.Ey,
|
||||
],
|
||||
]
|
||||
|
||||
for obj in slices:
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -673,16 +673,26 @@ class SubGridHSG(SubGridBaseGrid):
|
||||
|
||||
# Working region
|
||||
xs, ys, zs = self.round_to_grid(
|
||||
(self.i0 * self.dx * self.ratio, self.j0 * self.dy * self.ratio, self.k0 * self.dz * self.ratio)
|
||||
(
|
||||
self.i0 * self.dx * self.ratio,
|
||||
self.j0 * self.dy * self.ratio,
|
||||
self.k0 * self.dz * self.ratio,
|
||||
)
|
||||
)
|
||||
xf, yf, zf = self.round_to_grid(
|
||||
(self.i1 * self.dx * self.ratio, self.j1 * self.dy * self.ratio, self.k1 * self.dz * self.ratio)
|
||||
(
|
||||
self.i1 * self.dx * self.ratio,
|
||||
self.j1 * self.dy * self.ratio,
|
||||
self.k1 * self.dz * self.ratio,
|
||||
)
|
||||
)
|
||||
|
||||
logger.info("")
|
||||
logger.debug(f"[{self.name}] Type: {self.__class__.__name__}")
|
||||
logger.info(f"[{self.name}] Ratio: 1:{self.ratio}")
|
||||
logger.info(f"[{self.name}] Spatial discretisation: {self.dx:g} x " + f"{self.dy:g} x {self.dz:g}m")
|
||||
logger.info(
|
||||
f"[{self.name}] Spatial discretisation: {self.dx:g} x " + f"{self.dy:g} x {self.dz:g}m"
|
||||
)
|
||||
logger.info(
|
||||
f"[{self.name}] Extent (working region): {xs}m, {ys}m, {zs}m to {xf}m, {yf}m, {zf}m "
|
||||
+ f"(({self.nwx} x {self.nwy} x {self.nwz} = {self.nwx * self.nwy * self.nwz} cells)"
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Tobias Schruff
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -87,6 +87,27 @@ class CPUUpdates(Updates[GridType]):
|
||||
self.grid,
|
||||
)
|
||||
|
||||
# Update the magnetic and electric field components for the discrete plane wave
|
||||
for source in self.grid.discreteplanewaves:
|
||||
source.update_plane_wave(
|
||||
config.get_model_config().ompthreads,
|
||||
self.grid.updatecoeffsE,
|
||||
self.grid.updatecoeffsH,
|
||||
self.grid.Ex,
|
||||
self.grid.Ey,
|
||||
self.grid.Ez,
|
||||
self.grid.Hx,
|
||||
self.grid.Hy,
|
||||
self.grid.Hz,
|
||||
iteration,
|
||||
self.grid,
|
||||
cythonize=True,
|
||||
precompute=False,
|
||||
)
|
||||
|
||||
# if(self.grid.iteration % 10 == 0):
|
||||
# np.save('./snapshots/electric_z_{}.npy'.format(self.grid.iteration), self.grid.Ex)
|
||||
|
||||
def update_electric_a(self):
|
||||
"""Updates electric field components."""
|
||||
# All materials are non-dispersive so do standard update.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -124,7 +124,10 @@ class AddGrass(RotatableMixin, GeometryUserObject):
|
||||
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))
|
||||
fractalrange = (
|
||||
round_value(limits[0] / grid.dx),
|
||||
round_value(limits[1] / grid.dx),
|
||||
)
|
||||
# xminus surface
|
||||
if xs == volume.xs:
|
||||
logger.exception(
|
||||
@@ -151,7 +154,10 @@ class AddGrass(RotatableMixin, GeometryUserObject):
|
||||
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))
|
||||
fractalrange = (
|
||||
round_value(limits[0] / grid.dy),
|
||||
round_value(limits[1] / grid.dy),
|
||||
)
|
||||
# yminus surface
|
||||
if ys == volume.ys:
|
||||
logger.exception(
|
||||
@@ -175,7 +181,10 @@ class AddGrass(RotatableMixin, GeometryUserObject):
|
||||
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))
|
||||
fractalrange = (
|
||||
round_value(limits[0] / grid.dz),
|
||||
round_value(limits[1] / grid.dz),
|
||||
)
|
||||
# zminus surface
|
||||
if zs == volume.zs:
|
||||
logger.exception(
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -105,7 +105,7 @@ class AddSurfaceRoughness(RotatableMixin, GeometryUserObject):
|
||||
|
||||
if frac_dim < 0:
|
||||
logger.exception(
|
||||
f"{self.__str__()} requires a positive value for the " + "fractal dimension"
|
||||
f"{self.__str__()} requires a positive value for the fractal dimension"
|
||||
)
|
||||
raise ValueError
|
||||
if weighting[0] < 0:
|
||||
@@ -131,7 +131,10 @@ class AddSurfaceRoughness(RotatableMixin, GeometryUserObject):
|
||||
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))
|
||||
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:
|
||||
@@ -165,7 +168,10 @@ class AddSurfaceRoughness(RotatableMixin, GeometryUserObject):
|
||||
+ "surfaces of a fractal box"
|
||||
)
|
||||
raise ValueError
|
||||
fractalrange = (round_value(limits[0] / grid.dy), round_value(limits[1] / grid.dy))
|
||||
fractalrange = (
|
||||
round_value(limits[0] / grid.dy),
|
||||
round_value(limits[1] / grid.dy),
|
||||
)
|
||||
# yminus surface
|
||||
if ys == volume.ys:
|
||||
if fractalrange[0] < 0 or fractalrange[1] > volume.yf:
|
||||
@@ -196,7 +202,10 @@ class AddSurfaceRoughness(RotatableMixin, GeometryUserObject):
|
||||
+ "surfaces of a fractal box"
|
||||
)
|
||||
raise ValueError
|
||||
fractalrange = (round_value(limits[0] / grid.dz), round_value(limits[1] / grid.dz))
|
||||
fractalrange = (
|
||||
round_value(limits[0] / grid.dz),
|
||||
round_value(limits[1] / grid.dz),
|
||||
)
|
||||
# zminus surface
|
||||
if zs == volume.zs:
|
||||
if fractalrange[0] < 0 or fractalrange[1] > volume.zf:
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -147,7 +147,6 @@ class Box(RotatableMixin, GeometryUserObject):
|
||||
yf,
|
||||
zs,
|
||||
zf,
|
||||
config.get_model_config().ompthreads,
|
||||
numID,
|
||||
numIDx,
|
||||
numIDy,
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -665,7 +665,6 @@ class FractalBox(RotatableMixin, GeometryUserObject):
|
||||
self.volume.xs,
|
||||
self.volume.ys,
|
||||
self.volume.zs,
|
||||
config.get_model_config().ompthreads,
|
||||
waternumID,
|
||||
grassnumID,
|
||||
self.volume.averaging,
|
||||
@@ -700,7 +699,6 @@ class FractalBox(RotatableMixin, GeometryUserObject):
|
||||
self.volume.xs,
|
||||
self.volume.ys,
|
||||
self.volume.zs,
|
||||
config.get_model_config().ompthreads,
|
||||
0,
|
||||
self.volume.averaging,
|
||||
data,
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -151,7 +151,6 @@ class GeometryObjectsRead(GeometryUserObject):
|
||||
xs,
|
||||
ys,
|
||||
zs,
|
||||
config.get_model_config().ompthreads,
|
||||
numexistmaterials,
|
||||
averaging,
|
||||
data,
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -28,7 +28,6 @@ from scipy import interpolate
|
||||
|
||||
import gprMax.config as config
|
||||
from gprMax.grid.fdtd_grid import FDTDGrid
|
||||
from gprMax.grid.mpi_grid import MPIGrid
|
||||
from gprMax.materials import DispersiveMaterial as DispersiveMaterialUser
|
||||
from gprMax.materials import ListMaterial as ListMaterialUser
|
||||
from gprMax.materials import Material as MaterialUser
|
||||
@@ -36,13 +35,11 @@ from gprMax.materials import PeplinskiSoil as PeplinskiSoilUser
|
||||
from gprMax.materials import RangeMaterial as RangeMaterialUser
|
||||
from gprMax.pml import CFS, CFSParameter
|
||||
from gprMax.receivers import Rx as RxUser
|
||||
from gprMax.snapshots import MPISnapshot as MPISnapshotUser
|
||||
from gprMax.snapshots import Snapshot as SnapshotUser
|
||||
from gprMax.sources import DiscretePlaneWave as DiscretePlaneWaveUser
|
||||
from gprMax.sources import HertzianDipole as HertzianDipoleUser
|
||||
from gprMax.sources import MagneticDipole as MagneticDipoleUser
|
||||
from gprMax.sources import TransmissionLine as TransmissionLineUser
|
||||
from gprMax.sources import VoltageSource as VoltageSourceUser
|
||||
from gprMax.subgrids.grid import SubGridBaseGrid
|
||||
from gprMax.user_objects.cmds_geometry.cmds_geometry import (
|
||||
rotate_2point_object,
|
||||
rotate_polarisation,
|
||||
@@ -393,7 +390,7 @@ class VoltageSource(RotatableMixin, GridUserObject):
|
||||
x, y, z = uip.discretise_static_point(self.point)
|
||||
voltage_source.ID = f"{voltage_source.__class__.__name__}({x},{y},{z})"
|
||||
voltage_source.resistance = self.resistance
|
||||
voltage_source.waveform = grid.get_waveform_by_id(self.waveform_id)
|
||||
voltage_source.waveformID = self.waveform_id
|
||||
|
||||
if self.start is None or self.stop is None:
|
||||
voltage_source.start = 0
|
||||
@@ -402,7 +399,7 @@ class VoltageSource(RotatableMixin, GridUserObject):
|
||||
voltage_source.start = self.start
|
||||
voltage_source.stop = min(self.stop, grid.timewindow)
|
||||
|
||||
voltage_source.calculate_waveform_values(grid.iterations, grid.dt)
|
||||
voltage_source.calculate_waveform_values(grid)
|
||||
|
||||
return voltage_source
|
||||
|
||||
@@ -416,7 +413,7 @@ class VoltageSource(RotatableMixin, GridUserObject):
|
||||
f"{self.grid_name(grid)}Voltage source with polarity"
|
||||
f" {voltage_source.polarisation} at {x:g}m, {y:g}m, {z:g}m,"
|
||||
f" resistance {voltage_source.resistance:.1f} Ohms,"
|
||||
f"{startstop}using waveform {voltage_source.waveform.ID}"
|
||||
f"{startstop}using waveform {voltage_source.waveformID}"
|
||||
f" created."
|
||||
)
|
||||
|
||||
@@ -536,7 +533,7 @@ class HertzianDipole(RotatableMixin, GridUserObject):
|
||||
uip = self._create_uip(grid)
|
||||
x, y, z = uip.discretise_static_point(self.point)
|
||||
h.ID = f"{h.__class__.__name__}({x},{y},{z})"
|
||||
h.waveform = grid.get_waveform_by_id(self.waveform_id)
|
||||
h.waveformID = self.waveform_id
|
||||
|
||||
if self.start is None or self.stop is None:
|
||||
h.start = 0
|
||||
@@ -545,7 +542,7 @@ class HertzianDipole(RotatableMixin, GridUserObject):
|
||||
h.start = self.start
|
||||
h.stop = min(self.stop, grid.timewindow)
|
||||
|
||||
h.calculate_waveform_values(grid.iterations, grid.dt)
|
||||
h.calculate_waveform_values(grid)
|
||||
|
||||
return h
|
||||
|
||||
@@ -562,14 +559,14 @@ class HertzianDipole(RotatableMixin, GridUserObject):
|
||||
f"{self.grid_name(grid)}Hertzian dipole is a line source"
|
||||
f" in 2D with polarity {hertzian_dipole.polarisation}"
|
||||
f" at {x:g}m, {y:g}m, {z:g}m,{startstop}using"
|
||||
f" waveform {hertzian_dipole.waveform.ID} created."
|
||||
f" waveform {hertzian_dipole.waveformID} created."
|
||||
)
|
||||
else:
|
||||
logger.info(
|
||||
f"{self.grid_name(grid)}Hertzian dipole with polarity"
|
||||
f" {hertzian_dipole.polarisation} at {x:g}m, {y:g}m,"
|
||||
f" {z:g}m,{startstop}using waveform"
|
||||
f" {hertzian_dipole.waveform.ID} created."
|
||||
f" {hertzian_dipole.waveformID} created."
|
||||
)
|
||||
|
||||
def build(self, grid: FDTDGrid):
|
||||
@@ -694,7 +691,7 @@ class MagneticDipole(RotatableMixin, GridUserObject):
|
||||
uip = self._create_uip(grid)
|
||||
x, y, z = uip.discretise_static_point(self.point)
|
||||
m.ID = f"{m.__class__.__name__}({x},{y},{z})"
|
||||
m.waveform = grid.get_waveform_by_id(self.waveform_id)
|
||||
m.waveformID = self.waveform_id
|
||||
|
||||
if self.start is None or self.stop is None:
|
||||
m.start = 0
|
||||
@@ -703,7 +700,7 @@ class MagneticDipole(RotatableMixin, GridUserObject):
|
||||
m.start = self.start
|
||||
m.stop = min(self.stop, grid.timewindow)
|
||||
|
||||
m.calculate_waveform_values(grid.iterations, grid.dt)
|
||||
m.calculate_waveform_values(grid)
|
||||
|
||||
return m
|
||||
|
||||
@@ -716,7 +713,7 @@ class MagneticDipole(RotatableMixin, GridUserObject):
|
||||
logger.info(
|
||||
f"{self.grid_name(grid)}Magnetic dipole with polarity"
|
||||
f" {m.polarisation} at {x:g}m, {y:g}m, {z:g}m,"
|
||||
f"{startstop}using waveform {m.waveform.ID} created."
|
||||
f"{startstop}using waveform {m.waveformID} created."
|
||||
)
|
||||
|
||||
|
||||
@@ -850,7 +847,7 @@ class TransmissionLine(RotatableMixin, GridUserObject):
|
||||
x, y, z = uip.discretise_static_point(self.point)
|
||||
t.ID = f"{t.__class__.__name__}({x},{y},{z})"
|
||||
t.resistance = self.resistance
|
||||
t.waveform = grid.get_waveform_by_id(self.waveform_id)
|
||||
t.waveformID = self.waveform_id
|
||||
|
||||
if self.start is None or self.stop is None:
|
||||
t.start = 0
|
||||
@@ -859,7 +856,7 @@ class TransmissionLine(RotatableMixin, GridUserObject):
|
||||
t.start = self.start
|
||||
t.stop = min(self.stop, grid.timewindow)
|
||||
|
||||
t.calculate_waveform_values(grid.iterations, grid.dt)
|
||||
t.calculate_waveform_values(grid)
|
||||
t.calculate_incident_V_I(grid)
|
||||
|
||||
return t
|
||||
@@ -874,10 +871,160 @@ class TransmissionLine(RotatableMixin, GridUserObject):
|
||||
f"{self.grid_name(grid)}Transmission line with polarity"
|
||||
f" {t.polarisation} at {x:g}m, {y:g}m, {z:g}m,"
|
||||
f" resistance {t.resistance:.1f} Ohms,{startstop} using"
|
||||
f" waveform {t.waveform.ID} created."
|
||||
f" waveform {t.waveformID} created."
|
||||
)
|
||||
|
||||
|
||||
class DiscretePlaneWave(GridUserObject):
|
||||
"""
|
||||
Specifies a plane wave implemented using the discrete plane wave formulation.
|
||||
|
||||
Attributes:
|
||||
theta: float required for propagation angle (degrees) of wave.
|
||||
phi: float required for propagation angle (degrees) of wave.
|
||||
psi: float required for polarisation of wave.
|
||||
delta_theta: float optional for tolerance of theta angle to nearest
|
||||
rational angle.
|
||||
delta_phi: float optional for tolerance to phi angle to nearest
|
||||
rational angle.
|
||||
p1: tuple required for the lower left position (x, y, z) of the total
|
||||
field, scattered field (TFSF) box.
|
||||
p1: tuple required for the upper right position (x, y, z) of the total
|
||||
field, scattered field (TFSF) box.
|
||||
waveform_id: string required for identifier of waveform used with source.
|
||||
material_id: string optional of material identifier to use as the
|
||||
background material in the TFSF box.
|
||||
start: float optional to delay start time (secs) of source.
|
||||
stop: float optional to time (secs) to remove source.
|
||||
"""
|
||||
|
||||
@property
|
||||
def order(self):
|
||||
return 19
|
||||
|
||||
@property
|
||||
def hash(self):
|
||||
return "#discrete_plane_wave"
|
||||
|
||||
def __init__(self, **kwargs):
|
||||
super().__init__(**kwargs)
|
||||
|
||||
def build(self, grid: FDTDGrid):
|
||||
try:
|
||||
theta = self.kwargs["theta"]
|
||||
phi = self.kwargs["phi"]
|
||||
psi = self.kwargs["psi"]
|
||||
p1 = self.kwargs["p1"]
|
||||
p2 = self.kwargs["p2"]
|
||||
waveform_id = self.kwargs["waveform_id"]
|
||||
except KeyError:
|
||||
logger.exception(f"{self.params_str()} requires at least ten parameters.")
|
||||
raise
|
||||
try:
|
||||
dtheta = self.kwargs["delta_theta"]
|
||||
dphi = self.kwargs["delta_phi"]
|
||||
except KeyError:
|
||||
dtheta = 1.0
|
||||
dphi = 1.0
|
||||
|
||||
# Warn about using a discrete plane wave on GPU
|
||||
if config.sim_config.general["solver"] in ["cuda", "opencl"]:
|
||||
logger.exception(
|
||||
f"{self.params_str()} cannot currently be used "
|
||||
+ "with the CUDA or OpenCL-based solver. "
|
||||
)
|
||||
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(
|
||||
f"{self.params_str()} there is no waveform " + f"with the identifier {waveform_id}."
|
||||
)
|
||||
raise ValueError
|
||||
|
||||
if theta > 180:
|
||||
theta -= np.floor(theta / 180) * 180.0
|
||||
if phi > 360:
|
||||
phi -= np.floor(phi / 360) * 360.0
|
||||
if psi > 360:
|
||||
psi -= np.floor(psi / 360) * 360.0
|
||||
|
||||
uip = self._create_uip(grid)
|
||||
x_start, y_start, z_start = uip.check_src_rx_point(p1, self.params_str())
|
||||
x_stop, y_stop, z_stop = uip.check_src_rx_point(p2, self.params_str())
|
||||
|
||||
DPW = DiscretePlaneWaveUser(grid)
|
||||
DPW.corners = np.array([x_start, y_start, z_start, x_stop, y_stop, z_stop], dtype=np.int32)
|
||||
DPW.waveformID = waveform_id
|
||||
DPW.initializeDiscretePlaneWave(psi, phi, dphi, theta, dtheta, grid)
|
||||
|
||||
try:
|
||||
DPW.material_ID = self.kwargs["material_id"]
|
||||
except KeyError:
|
||||
DPW.material_ID = 1
|
||||
|
||||
try:
|
||||
# Check source start & source remove time parameters
|
||||
start = self.kwargs["start"]
|
||||
stop = self.kwargs["stop"]
|
||||
if start < 0:
|
||||
logger.exception(
|
||||
self.params_str()
|
||||
+ (" 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.")
|
||||
)
|
||||
raise ValueError
|
||||
if stop - start <= 0:
|
||||
logger.exception(
|
||||
self.params_str() + (" duration of the source should not be zero or less.")
|
||||
)
|
||||
raise ValueError
|
||||
DPW.start = start
|
||||
DPW.stop = min(stop, grid.timewindow)
|
||||
startstop = f" start time {t.start:g} secs, finish time {t.stop:g} secs "
|
||||
except KeyError:
|
||||
DPW.start = 0
|
||||
DPW.stop = grid.timewindow
|
||||
startstop = " "
|
||||
|
||||
precompute = True
|
||||
if precompute:
|
||||
DPW.calculate_waveform_values(grid)
|
||||
|
||||
logger.info(
|
||||
f"{self.grid_name(grid)}Discrete Plane Wave within the TFSF Box "
|
||||
+ f"spanning from {p1} m to {p2} m, incident in the direction "
|
||||
+ f"theta {theta} degrees and phi {phi} degrees "
|
||||
+ startstop
|
||||
+ f"using waveform {DPW.waveformID} created."
|
||||
)
|
||||
phi_approx = np.arctan2(DPW.m[1] / grid.dy, DPW.m[0] / grid.dx) * 180 / np.pi
|
||||
theta_approx = (
|
||||
np.arctan2(
|
||||
np.sqrt(
|
||||
(DPW.m[0] / grid.dx) * (DPW.m[0] / grid.dx)
|
||||
+ (DPW.m[1] / grid.dy) * (DPW.m[1] / grid.dy)
|
||||
),
|
||||
DPW.m[2] / grid.dz,
|
||||
)
|
||||
* 180
|
||||
/ np.pi
|
||||
)
|
||||
logger.info(
|
||||
f"{self.grid_name(grid)}Discrete Plane Wave has been discretized "
|
||||
+ "the angles have been approximated to the nearest rational angles "
|
||||
+ "with some small tolerance levels. The chosen rational integers are "
|
||||
+ f"[m_x, m_y, m_z] : {DPW.m[:3]}. The approximated angles are: "
|
||||
+ f"Phi: {phi_approx:.3f} and Theta: {theta_approx:.3f}"
|
||||
)
|
||||
|
||||
grid.discreteplanewaves.append(DPW)
|
||||
|
||||
|
||||
class Rx(RotatableMixin, GridUserObject):
|
||||
"""Specifies output points in the model.
|
||||
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -170,7 +170,9 @@ def get_host_info():
|
||||
)
|
||||
model = (
|
||||
subprocess.check_output(
|
||||
["cat", "/sys/class/dmi/id/product_name"], shell=False, stderr=subprocess.STDOUT
|
||||
["cat", "/sys/class/dmi/id/product_name"],
|
||||
shell=False,
|
||||
stderr=subprocess.STDOUT,
|
||||
)
|
||||
.decode("utf-8")
|
||||
.strip()
|
||||
@@ -248,7 +250,7 @@ def print_host_info(hostinfo):
|
||||
"""
|
||||
|
||||
hyperthreadingstr = (
|
||||
f", {config.sim_config.hostinfo['logicalcores']} " f"cores with Hyper-Threading"
|
||||
f", {config.sim_config.hostinfo['logicalcores']} cores with Hyper-Threading"
|
||||
if config.sim_config.hostinfo["hyperthreading"]
|
||||
else ""
|
||||
)
|
||||
@@ -506,12 +508,12 @@ def detect_cuda_gpus():
|
||||
|
||||
|
||||
def print_cuda_info(devs):
|
||||
""""Prints info about detected CUDA-capable GPU(s).
|
||||
"""Prints info about detected CUDA-capable GPU(s).
|
||||
|
||||
Args:
|
||||
devs: dict of detected pycuda device object(s) where where device ID(s)
|
||||
are keys.
|
||||
""" ""
|
||||
"""
|
||||
|
||||
import pycuda
|
||||
|
||||
@@ -561,12 +563,12 @@ def detect_opencl():
|
||||
|
||||
|
||||
def print_opencl_info(devs):
|
||||
""""Prints info about detected OpenCL-capable device(s).
|
||||
"""Prints info about detected OpenCL-capable device(s).
|
||||
|
||||
Args:
|
||||
devs: dict of detected pyopencl device object(s) where where device ID(s)
|
||||
are keys.
|
||||
""" ""
|
||||
"""
|
||||
|
||||
import pyopencl as cl
|
||||
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -22,6 +22,7 @@ import logging
|
||||
import re
|
||||
import textwrap
|
||||
from shutil import get_terminal_size
|
||||
from time import perf_counter as timer_fn
|
||||
from typing import Union
|
||||
|
||||
import numpy as np
|
||||
@@ -31,12 +32,6 @@ init()
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
try:
|
||||
from time import thread_time as timer_fn
|
||||
except ImportError:
|
||||
# thread_time not always available in macOS
|
||||
from time import process_time as timer_fn
|
||||
|
||||
|
||||
def get_terminal_width():
|
||||
"""Gets/sets width of terminal being used.
|
||||
@@ -100,7 +95,7 @@ def logo(version):
|
||||
)
|
||||
|
||||
logo = (
|
||||
""" www.gprmax.com __ __
|
||||
r""" www.gprmax.com __ __
|
||||
__ _ _ __ _ __| \/ | __ ___ __
|
||||
/ _` | '_ \| '__| |\/| |/ _` \ \/ /
|
||||
| (_| | |_) | | | | | | (_| |> <
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -66,7 +66,13 @@ class Waveform:
|
||||
waveforms.
|
||||
"""
|
||||
|
||||
if self.type in ["gaussian", "gaussiandot", "gaussiandotnorm", "gaussianprime", "gaussiandoubleprime"]:
|
||||
if self.type in [
|
||||
"gaussian",
|
||||
"gaussiandot",
|
||||
"gaussiandotnorm",
|
||||
"gaussianprime",
|
||||
"gaussiandoubleprime",
|
||||
]:
|
||||
self.chi = 1 / self.freq
|
||||
self.zeta = 2 * np.pi**2 * self.freq**2
|
||||
elif self.type in ["gaussiandotdot", "gaussiandotdotnorm", "ricker"]:
|
||||
@@ -102,18 +108,27 @@ class Waveform:
|
||||
|
||||
elif self.type in ["gaussiandotdot", "gaussiandoubleprime"]:
|
||||
delay = time - self.chi
|
||||
ampvalue = 2 * self.zeta * (2 * self.zeta * delay**2 - 1) * np.exp(-self.zeta * delay**2)
|
||||
ampvalue = (
|
||||
2 * self.zeta * (2 * self.zeta * delay**2 - 1) * np.exp(-self.zeta * delay**2)
|
||||
)
|
||||
|
||||
elif self.type == "gaussiandotdotnorm":
|
||||
delay = time - self.chi
|
||||
normalise = 1 / (2 * self.zeta)
|
||||
ampvalue = 2 * self.zeta * (2 * self.zeta * delay**2 - 1) * np.exp(-self.zeta * delay**2) * normalise
|
||||
ampvalue = (
|
||||
2
|
||||
* self.zeta
|
||||
* (2 * self.zeta * delay**2 - 1)
|
||||
* np.exp(-self.zeta * delay**2)
|
||||
* normalise
|
||||
)
|
||||
|
||||
elif self.type == "ricker":
|
||||
delay = time - self.chi
|
||||
normalise = 1 / (2 * self.zeta)
|
||||
ampvalue = -(
|
||||
(2 * self.zeta * (2 * self.zeta * delay**2 - 1) * np.exp(-self.zeta * delay**2)) * normalise
|
||||
(2 * self.zeta * (2 * self.zeta * delay**2 - 1) * np.exp(-self.zeta * delay**2))
|
||||
* normalise
|
||||
)
|
||||
|
||||
elif self.type == "sine":
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
[build-system]
|
||||
requires = ['setuptools', 'wheel', 'numpy>=1.19.0', 'Cython>=0.29.21',
|
||||
requires = ['setuptools>=64', 'wheel', 'numpy>=1.19.0', 'Cython>=0.29.21',
|
||||
'jinja2']
|
||||
build-backend = "setuptools.build_meta"
|
||||
|
||||
二进制文件未显示。
二进制文件未显示。
@@ -0,0 +1,8 @@
|
||||
#title: Discrete plane wave in free-space
|
||||
#domain: 0.100 0.100 0.100
|
||||
#dx_dy_dz: 0.001 0.001 0.001
|
||||
#time_window: 3e-9
|
||||
|
||||
#waveform: gaussiandot 1 1e9 myWave
|
||||
#discrete_plane_wave: 0.010 0.010 0.010 0.040 0.040 0.040 90.0 63.4 36.7 myWave
|
||||
#rx: 0.070 0.070 0.070
|
||||
@@ -22,6 +22,13 @@ class TestMagneticDipoleSource(GprMaxRegressionTest):
|
||||
model = parameter(["magnetic_dipole_fs"])
|
||||
|
||||
|
||||
@rfm.simple_test
|
||||
class TestDiscretePlaneWaveSource(GprMaxRegressionTest):
|
||||
tags = {"test", "serial", "discrete_plane_wave", "waveform"}
|
||||
sourcesdir = "src/source_tests"
|
||||
model = parameter(["discrete_plane_wave_fs"])
|
||||
|
||||
|
||||
@rfm.simple_test
|
||||
class TestTransmissionLineSource(GprMaxRegressionTest):
|
||||
tags = {"test", "serial", "transmission_line", "waveform"}
|
||||
|
||||
29
setup.py
29
setup.py
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -33,7 +33,9 @@ from setuptools import Extension, find_packages, setup
|
||||
# Check Python version
|
||||
MIN_PYTHON_VERSION = (3, 7)
|
||||
if sys.version_info[:2] < MIN_PYTHON_VERSION:
|
||||
sys.exit("\nExited: Requires Python {MIN_PYTHON_VERSION[0]}.{MIN_PYTHON_VERSION[1]} or newer!\n")
|
||||
sys.exit(
|
||||
"\nExited: Requires Python {MIN_PYTHON_VERSION[0]}.{MIN_PYTHON_VERSION[1]} or newer!\n"
|
||||
)
|
||||
|
||||
# Importing gprMax _version__.py before building can cause issues.
|
||||
with open("gprMax/_version.py", "r") as fd:
|
||||
@@ -144,9 +146,9 @@ if "cleanall" in sys.argv:
|
||||
except OSError:
|
||||
print(f"Could not remove: {filebase}.c")
|
||||
# Remove compiled Cython modules
|
||||
libfile = glob.glob(os.path.join(os.getcwd(), os.path.splitext(file)[0]) + "*.pyd") + glob.glob(
|
||||
os.path.join(os.getcwd(), os.path.splitext(file)[0]) + "*.so"
|
||||
)
|
||||
libfile = glob.glob(
|
||||
os.path.join(os.getcwd(), os.path.splitext(file)[0]) + "*.pyd"
|
||||
) + glob.glob(os.path.join(os.getcwd(), os.path.splitext(file)[0]) + "*.so")
|
||||
if libfile:
|
||||
libfile = libfile[0]
|
||||
try:
|
||||
@@ -182,7 +184,11 @@ else:
|
||||
elif sys.platform == "darwin":
|
||||
# Check for Intel or Apple M series CPU
|
||||
cpuID = (
|
||||
subprocess.check_output("sysctl -n machdep.cpu.brand_string", shell=True, stderr=subprocess.STDOUT)
|
||||
subprocess.check_output(
|
||||
"sysctl -n machdep.cpu.brand_string",
|
||||
shell=True,
|
||||
stderr=subprocess.STDOUT,
|
||||
)
|
||||
.decode("utf-8")
|
||||
.strip()
|
||||
)
|
||||
@@ -214,7 +220,13 @@ else:
|
||||
pass
|
||||
os.environ["MIN_SUPPORTED_MACOSX_DEPLOYMENT_TARGET"] = MIN_MACOS_VERSION
|
||||
# Sometimes worth testing with '-fstrict-aliasing', '-fno-common'
|
||||
compile_args = ["-O3", "-w", "-fopenmp", "-march=native", f"-mmacosx-version-min={MIN_MACOS_VERSION}"]
|
||||
compile_args = [
|
||||
"-O3",
|
||||
"-w",
|
||||
"-fopenmp",
|
||||
"-march=native",
|
||||
f"-mmacosx-version-min={MIN_MACOS_VERSION}",
|
||||
]
|
||||
linker_args = ["-fopenmp", f"-mmacosx-version-min={MIN_MACOS_VERSION}"]
|
||||
libraries = ["gomp"]
|
||||
|
||||
@@ -261,7 +273,8 @@ else:
|
||||
version=version,
|
||||
author="Craig Warren, Antonis Giannopoulos, and John Hartley",
|
||||
url="http://www.gprmax.com",
|
||||
description="Electromagnetic Modelling Software based on the " + "Finite-Difference Time-Domain (FDTD) method",
|
||||
description="Electromagnetic Modelling Software based on the "
|
||||
+ "Finite-Difference Time-Domain (FDTD) method",
|
||||
long_description=long_description,
|
||||
long_description_content_type="text/x-rst",
|
||||
license="GPLv3+",
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -164,7 +164,9 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx):
|
||||
tmp = Hy_x / Hy_y
|
||||
except ZeroDivisionError:
|
||||
tmp = 0
|
||||
fields[timestep, 4] = -tmp * (-(Hy_y / (4 * np.pi * Hr_y**3)) * (f_Hy + (tau_Hy * fdot_Hy)))
|
||||
fields[timestep, 4] = -tmp * (
|
||||
-(Hy_y / (4 * np.pi * Hr_y**3)) * (f_Hy + (tau_Hy * fdot_Hy))
|
||||
)
|
||||
|
||||
# Hz
|
||||
fields[timestep, 5] = 0
|
||||
|
||||
@@ -1,9 +1,8 @@
|
||||
"""A series of models with different domain sizes used for benchmarking.
|
||||
The domain is free space with a simple source (Hertzian Dipole) and
|
||||
receiver at the centre.
|
||||
The domain is free space with a simple source (Hertzian Dipole) and
|
||||
receiver at the centre.
|
||||
"""
|
||||
|
||||
|
||||
import itertools
|
||||
from pathlib import Path
|
||||
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -82,7 +82,10 @@ def diff_output_files(filename1, filename2):
|
||||
for i in range(len(outputs2)):
|
||||
maxi = np.amax(np.abs(data1[:, i]))
|
||||
datadiffs[:, i] = np.divide(
|
||||
np.abs(data2[:, i] - data1[:, i]), maxi, out=np.zeros_like(data1[:, i]), where=maxi != 0
|
||||
np.abs(data2[:, i] - data1[:, i]),
|
||||
maxi,
|
||||
out=np.zeros_like(data1[:, i]),
|
||||
where=maxi != 0,
|
||||
) # Replace any division by zero with zero
|
||||
|
||||
# Calculate power (ignore warning from taking a log of any zero values)
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -46,12 +46,16 @@ maxerrors = []
|
||||
testmodels = [basename + "_" + s for s in PMLIDs]
|
||||
|
||||
fig, ax = plt.subplots(
|
||||
subplot_kw=dict(xlabel="Iterations", ylabel="Error [dB]"), figsize=(20, 10), facecolor="w", edgecolor="w"
|
||||
subplot_kw=dict(xlabel="Iterations", ylabel="Error [dB]"),
|
||||
figsize=(20, 10),
|
||||
facecolor="w",
|
||||
edgecolor="w",
|
||||
)
|
||||
|
||||
for x, model in enumerate(testmodels):
|
||||
time, datadiffs = diff_output_files(
|
||||
fn.parent.joinpath(basename + "_ref.h5"), fn.parent.joinpath(basename + str(x + 1) + ".h5")
|
||||
fn.parent.joinpath(basename + "_ref.h5"),
|
||||
fn.parent.joinpath(basename + str(x + 1) + ".h5"),
|
||||
)
|
||||
|
||||
# Print maximum error value
|
||||
@@ -60,7 +64,14 @@ for x, model in enumerate(testmodels):
|
||||
logger.info(f"{model}: Max. error {maxerrors[x]}")
|
||||
|
||||
# Plot diffs (select column to choose field component, 0-Ex, 1-Ey etc..)
|
||||
ax.plot(time[start::], datadiffs[start::, 1], color=next(colors), lw=2, ls=next(lines), label=model)
|
||||
ax.plot(
|
||||
time[start::],
|
||||
datadiffs[start::, 1],
|
||||
color=next(colors),
|
||||
lw=2,
|
||||
ls=next(lines),
|
||||
label=model,
|
||||
)
|
||||
ax.set_xticks(np.arange(0, 2200, step=100))
|
||||
ax.set_xlim([0, 2100])
|
||||
ax.set_yticks(np.arange(-160, 0, step=20))
|
||||
|
||||
@@ -22,7 +22,9 @@ time_window = gprMax.TimeWindow(iterations=2100)
|
||||
tssf = gprMax.TimeStepStabilityFactor(f=0.99)
|
||||
|
||||
waveform = gprMax.Waveform(wave_type="gaussiandotnorm", amp=1, freq=9.42e9, id="mypulse")
|
||||
hertzian_dipole = gprMax.HertzianDipole(polarisation="z", p1=(0.013, 0.013, 0.014), waveform_id="mypulse")
|
||||
hertzian_dipole = gprMax.HertzianDipole(
|
||||
polarisation="z", p1=(0.013, 0.013, 0.014), waveform_id="mypulse"
|
||||
)
|
||||
rx = gprMax.Rx(p1=(0.038, 0.114, 0.013))
|
||||
|
||||
plate = gprMax.Plate(p1=(0.013, 0.013, 0.013), p2=(0.038, 0.113, 0.013), material_id="pec")
|
||||
|
||||
@@ -22,7 +22,9 @@ time_window = gprMax.TimeWindow(iterations=2100)
|
||||
tssf = gprMax.TimeStepStabilityFactor(f=0.99)
|
||||
|
||||
waveform = gprMax.Waveform(wave_type="gaussiandotnorm", amp=1, freq=9.42e9, id="mypulse")
|
||||
hertzian_dipole = gprMax.HertzianDipole(polarisation="z", p1=(0.088, 0.088, 0.089), waveform_id="mypulse")
|
||||
hertzian_dipole = gprMax.HertzianDipole(
|
||||
polarisation="z", p1=(0.088, 0.088, 0.089), waveform_id="mypulse"
|
||||
)
|
||||
rx = gprMax.Rx(p1=(0.113, 0.189, 0.088))
|
||||
|
||||
plate = gprMax.Plate(p1=(0.088, 0.088, 0.088), p2=(0.113, 0.188, 0.088), material_id="pec")
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -74,7 +74,13 @@ for x, PMLID in enumerate(PMLIDs):
|
||||
ax.set_ylabel(f"{output} error [dB]")
|
||||
|
||||
# Save a PDF/PNG of the figure
|
||||
fig.savefig(basename + "_diffs_" + PMLID + ".pdf", dpi=None, format="pdf", bbox_inches="tight", pad_inches=0.1)
|
||||
fig.savefig(
|
||||
basename + "_diffs_" + PMLID + ".pdf",
|
||||
dpi=None,
|
||||
format="pdf",
|
||||
bbox_inches="tight",
|
||||
pad_inches=0.1,
|
||||
)
|
||||
# fig.savefig(basename + "_diffs_" + PMLID + ".png", dpi=150, format='png', bbox_inches='tight', pad_inches=0.1)
|
||||
|
||||
plt.show()
|
||||
|
||||
@@ -18,7 +18,9 @@ dxdydz = gprMax.Discretisation(p1=(dl, dl, dl))
|
||||
time_window = gprMax.TimeWindow(time=3e-9)
|
||||
|
||||
waveform = gprMax.Waveform(wave_type="gaussian", amp=1, freq=1e9, id="mypulse")
|
||||
hertzian_dipole = gprMax.HertzianDipole(polarisation="z", p1=(0.050, 0.050, 0.050), waveform_id="mypulse")
|
||||
hertzian_dipole = gprMax.HertzianDipole(
|
||||
polarisation="z", p1=(0.050, 0.050, 0.050), waveform_id="mypulse"
|
||||
)
|
||||
rx = gprMax.Rx(p1=(0.070, 0.070, 0.070))
|
||||
|
||||
# PML cases
|
||||
|
||||
@@ -31,7 +31,10 @@ tx_pos = (x / 2, y / 2, z / 2)
|
||||
# Source excitation and type
|
||||
wave = gprMax.Waveform(wave_type="gaussian", amp=1, freq=1.5e9, id="mypulse")
|
||||
tl = gprMax.TransmissionLine(
|
||||
p1=(tx_pos[0], tx_pos[1], tx_pos[2]), polarisation="x", resistance=50, waveform_id="mypulse"
|
||||
p1=(tx_pos[0], tx_pos[1], tx_pos[2]),
|
||||
polarisation="x",
|
||||
resistance=50,
|
||||
waveform_id="mypulse",
|
||||
)
|
||||
scene.add(wave)
|
||||
scene.add(tl)
|
||||
@@ -59,8 +62,16 @@ scene.add(t2)
|
||||
|
||||
# Detailed geometry view around bowtie and feed position
|
||||
gv1 = gprMax.GeometryView(
|
||||
p1=(tx_pos[0] - bowtie_dims[0] - 2 * dl, tx_pos[1] - bowtie_dims[1] / 2 - 2 * dl, tx_pos[2] - 2 * dl),
|
||||
p2=(tx_pos[0] + bowtie_dims[0] + 2 * dl, tx_pos[1] + bowtie_dims[1] / 2 + 2 * dl, tx_pos[2] + 2 * dl),
|
||||
p1=(
|
||||
tx_pos[0] - bowtie_dims[0] - 2 * dl,
|
||||
tx_pos[1] - bowtie_dims[1] / 2 - 2 * dl,
|
||||
tx_pos[2] - 2 * dl,
|
||||
),
|
||||
p2=(
|
||||
tx_pos[0] + bowtie_dims[0] + 2 * dl,
|
||||
tx_pos[1] + bowtie_dims[1] / 2 + 2 * dl,
|
||||
tx_pos[2] + 2 * dl,
|
||||
),
|
||||
dl=(dl, dl, dl),
|
||||
filename="antenna_bowtie_fs_pcb",
|
||||
output_type="f",
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -32,7 +32,8 @@ logger = logging.getLogger(__name__)
|
||||
|
||||
# Parse command line arguments
|
||||
parser = argparse.ArgumentParser(
|
||||
description="Plots a comparison of fields between " + "given simulation output and experimental data files.",
|
||||
description="Plots a comparison of fields between "
|
||||
+ "given simulation output and experimental data files.",
|
||||
usage="cd gprMax; python -m testing.test_experimental modelfile realfile output",
|
||||
)
|
||||
parser.add_argument("modelfile", help="name of model output file including path")
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# Copyright (C) 2015-2024: The University of Edinburgh, United Kingdom
|
||||
# Copyright (C) 2015-2025: The University of Edinburgh, United Kingdom
|
||||
# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley
|
||||
#
|
||||
# This file is part of gprMax.
|
||||
@@ -113,7 +113,11 @@ for i, model in enumerate(testmodels):
|
||||
# Tx/Rx position to feed to analytical solution
|
||||
rxpos = filetest[path].attrs["Position"]
|
||||
txpos = filetest["/srcs/src1/"].attrs["Position"]
|
||||
rxposrelative = ((rxpos[0] - txpos[0]), (rxpos[1] - txpos[1]), (rxpos[2] - txpos[2]))
|
||||
rxposrelative = (
|
||||
(rxpos[0] - txpos[0]),
|
||||
(rxpos[1] - txpos[1]),
|
||||
(rxpos[2] - txpos[2]),
|
||||
)
|
||||
|
||||
# Analytical solution of a dipole in free space
|
||||
dataref = hertzian_dipole_fs(
|
||||
|
||||
@@ -11,12 +11,10 @@ import os
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
|
||||
from scipy.constants import c
|
||||
from scipy.constants import epsilon_0 as e0
|
||||
from scipy.constants import mu_0 as m0
|
||||
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
|
||||
@@ -140,7 +138,11 @@ ax.set_yticklabels(yticks)
|
||||
ax.grid(True)
|
||||
handles, existlabels = ax.get_legend_handles_labels()
|
||||
leg = ax.legend(
|
||||
[handles[0], handles[-1]], [existlabels[0], existlabels[-1]], ncol=2, loc=(0.27, -0.12), frameon=False
|
||||
[handles[0], handles[-1]],
|
||||
[existlabels[0], existlabels[-1]],
|
||||
ncol=2,
|
||||
loc=(0.27, -0.12),
|
||||
frameon=False,
|
||||
) # Plot just first and last legend entries
|
||||
# leg = ax.legend([handles[0], handles[-3], handles[-2], handles[-1]], [existlabels[0], existlabels[-3], existlabels[-2], existlabels[-1]], ncol=4, loc=(-0.13,-0.12), frameon=False)
|
||||
[legobj.set_linewidth(2) for legobj in leg.legendHandles]
|
||||
|
||||
@@ -9,9 +9,13 @@ logger = logging.getLogger(__name__)
|
||||
|
||||
# Parse command line arguments
|
||||
parser = argparse.ArgumentParser(
|
||||
description="Writes a HDF5 file of AustinMan or AustinWoman head only.", usage="python head_only_hdf5 filename"
|
||||
description="Writes a HDF5 file of AustinMan or AustinWoman head only.",
|
||||
usage="python head_only_hdf5 filename",
|
||||
)
|
||||
parser.add_argument(
|
||||
"filename",
|
||||
help="name and path to (HDF5) file containing AustinMan or AustinWoman model",
|
||||
)
|
||||
parser.add_argument("filename", help="name and path to (HDF5) file containing AustinMan or AustinWoman model")
|
||||
args = parser.parse_args()
|
||||
|
||||
# Read full body HDF5 file
|
||||
@@ -22,7 +26,9 @@ data = f["/data"][:, :, :]
|
||||
# Define head as last 1/8 of total body height
|
||||
nzhead = 7 * int(data.shape[2] / 8)
|
||||
|
||||
logger.info(f"Dimensions of head model: {data.shape[0]:g} x {data.shape[1]:g} x {data.shape[2] - nzhead:g} cells")
|
||||
logger.info(
|
||||
f"Dimensions of head model: {data.shape[0]:g} x {data.shape[1]:g} x {data.shape[2] - nzhead:g} cells"
|
||||
)
|
||||
|
||||
# Write HDF5 file
|
||||
headfile = os.path.splitext(args.filename)[0] + "_head.h5"
|
||||
|
||||
@@ -100,7 +100,15 @@ class Relaxation(object):
|
||||
def check_inputs(self):
|
||||
"""Check the validity of the inputs."""
|
||||
try:
|
||||
d = [float(i) for i in [self.number_of_debye_poles, self.sigma, self.mu, self.mu_sigma]]
|
||||
d = [
|
||||
float(i)
|
||||
for i in [
|
||||
self.number_of_debye_poles,
|
||||
self.sigma,
|
||||
self.mu,
|
||||
self.mu_sigma,
|
||||
]
|
||||
]
|
||||
except ValueError:
|
||||
sys.exit("The inputs should be numeric.")
|
||||
if not isinstance(self.number_of_debye_poles, int):
|
||||
@@ -120,7 +128,7 @@ class Relaxation(object):
|
||||
Returns:
|
||||
s (str): Info about chosen function and its parameters.
|
||||
"""
|
||||
print(f"Approximating {self.name}" f" using {self.number_of_debye_poles} Debye poles")
|
||||
print(f"Approximating {self.name} using {self.number_of_debye_poles} Debye poles")
|
||||
print(f"{self.name} parameters: ")
|
||||
s = "".join(f"{k:10s} = {v}\n" for k, v in self.params.items())
|
||||
print(s)
|
||||
@@ -172,7 +180,12 @@ class Relaxation(object):
|
||||
self.rl, self.im = q.real, q.imag
|
||||
|
||||
if self.number_of_debye_poles == -1:
|
||||
print("\n#########", "Try to automaticaly fit number of Debye poles, up to 20!", "##########\n", sep="")
|
||||
print(
|
||||
"\n#########",
|
||||
"Try to automaticaly fit number of Debye poles, up to 20!",
|
||||
"##########\n",
|
||||
sep="",
|
||||
)
|
||||
error = np.infty # artificial best error starting value
|
||||
self.number_of_debye_poles = 1
|
||||
iteration = 1
|
||||
@@ -195,7 +208,11 @@ class Relaxation(object):
|
||||
|
||||
# Print the results in gprMax format style
|
||||
properties = self.print_output(tau, weights, ee)
|
||||
print(f"The average fractional error for:\n" f"- real part: {err_real}\n" f"- imaginary part: {err_imag}\n")
|
||||
print(
|
||||
f"The average fractional error for:\n"
|
||||
f"- real part: {err_real}\n"
|
||||
f"- imaginary part: {err_imag}\n"
|
||||
)
|
||||
if self.save:
|
||||
self.save_result(properties)
|
||||
# Plot the actual and the approximate dielectric properties
|
||||
@@ -225,11 +242,13 @@ class Relaxation(object):
|
||||
|
||||
# Print the Debye expnasion in a gprMax format
|
||||
material_prop = []
|
||||
material_prop.append(f"#material: {ee} {self.sigma} {self.mu} {self.mu_sigma} {self.material_name}\n")
|
||||
material_prop.append(
|
||||
f"#material: {ee} {self.sigma} {self.mu} {self.mu_sigma} {self.material_name}\n"
|
||||
)
|
||||
print(material_prop[0], end="")
|
||||
dispersion_prop = f"#add_dispersion_debye: {len(tau)}"
|
||||
for i in range(len(tau)):
|
||||
dispersion_prop += f" {weights[i]} {10**tau[i]}"
|
||||
dispersion_prop += f" {weights[i]} {10 ** tau[i]}"
|
||||
dispersion_prop += f" {self.material_name}"
|
||||
print(dispersion_prop)
|
||||
material_prop.append(dispersion_prop + "\n")
|
||||
@@ -251,10 +270,34 @@ class Relaxation(object):
|
||||
gs = gridspec.GridSpec(2, 1)
|
||||
ax = fig.add_subplot(gs[0])
|
||||
ax.grid(b=True, which="major", linewidth=0.2, linestyle="--")
|
||||
ax.semilogx(self.freq * 1e-6, rl_exp, "b-", linewidth=2.0, label="Debye Expansion: Real part")
|
||||
ax.semilogx(self.freq * 1e-6, -im_exp, "k-", linewidth=2.0, label="Debye Expansion: Imaginary part")
|
||||
ax.semilogx(self.freq * 1e-6, self.rl, "r.", linewidth=2.0, label=f"{self.name}: Real part")
|
||||
ax.semilogx(self.freq * 1e-6, -self.im, "g.", linewidth=2.0, label=f"{self.name}: Imaginary part")
|
||||
ax.semilogx(
|
||||
self.freq * 1e-6,
|
||||
rl_exp,
|
||||
"b-",
|
||||
linewidth=2.0,
|
||||
label="Debye Expansion: Real part",
|
||||
)
|
||||
ax.semilogx(
|
||||
self.freq * 1e-6,
|
||||
-im_exp,
|
||||
"k-",
|
||||
linewidth=2.0,
|
||||
label="Debye Expansion: Imaginary part",
|
||||
)
|
||||
ax.semilogx(
|
||||
self.freq * 1e-6,
|
||||
self.rl,
|
||||
"r.",
|
||||
linewidth=2.0,
|
||||
label=f"{self.name}: Real part",
|
||||
)
|
||||
ax.semilogx(
|
||||
self.freq * 1e-6,
|
||||
-self.im,
|
||||
"g.",
|
||||
linewidth=2.0,
|
||||
label=f"{self.name}: Imaginary part",
|
||||
)
|
||||
ax.set_ylim([-1, np.max(np.concatenate([self.rl, -self.im])) + 1])
|
||||
ax.legend()
|
||||
ax.set_xlabel("Frequency (MHz)")
|
||||
@@ -262,8 +305,20 @@ class Relaxation(object):
|
||||
|
||||
ax = fig.add_subplot(gs[1])
|
||||
ax.grid(b=True, which="major", linewidth=0.2, linestyle="--")
|
||||
ax.semilogx(self.freq * 1e-6, (rl_exp - self.rl) / (self.rl + 1), "b-", linewidth=2.0, label="Real part")
|
||||
ax.semilogx(self.freq * 1e-6, (-im_exp + self.im) / (self.im + 1), "k-", linewidth=2.0, label="Imaginary part")
|
||||
ax.semilogx(
|
||||
self.freq * 1e-6,
|
||||
(rl_exp - self.rl) / (self.rl + 1),
|
||||
"b-",
|
||||
linewidth=2.0,
|
||||
label="Real part",
|
||||
)
|
||||
ax.semilogx(
|
||||
self.freq * 1e-6,
|
||||
(-im_exp + self.im) / (self.im + 1),
|
||||
"k-",
|
||||
linewidth=2.0,
|
||||
label="Imaginary part",
|
||||
)
|
||||
ax.legend()
|
||||
ax.set_xlabel("Frequency (MHz)")
|
||||
ax.set_ylabel("Relative approximation error")
|
||||
@@ -306,7 +361,9 @@ class Relaxation(object):
|
||||
elif os.path.isdir("user_libs/materials"):
|
||||
file_path = os.path.join("user_libs", "materials", "my_materials.txt")
|
||||
else:
|
||||
sys.exit("Cannot save material properties " f"in {os.path.join(fdir, 'my_materials.txt')}!")
|
||||
sys.exit(
|
||||
"Cannot save material properties " f"in {os.path.join(fdir, 'my_materials.txt')}!"
|
||||
)
|
||||
with open(file_path, "a") as fileH:
|
||||
fileH.write(f"## {output[0].split(' ')[-1]}")
|
||||
fileH.writelines(output)
|
||||
@@ -382,7 +439,13 @@ class HavriliakNegami(Relaxation):
|
||||
self.f_min, self.f_max = f_min, f_max
|
||||
# Choosing n frequencies logarithmicaly equally spaced between the bounds given
|
||||
self.set_freq(self.f_min, self.f_max, self.f_n)
|
||||
self.e_inf, self.alpha, self.beta, self.de, self.tau_0 = e_inf, alpha, beta, de, tau_0
|
||||
self.e_inf, self.alpha, self.beta, self.de, self.tau_0 = (
|
||||
e_inf,
|
||||
alpha,
|
||||
beta,
|
||||
de,
|
||||
tau_0,
|
||||
)
|
||||
self.params = {
|
||||
"f_min": self.f_min,
|
||||
"f_max": self.f_max,
|
||||
@@ -412,7 +475,10 @@ class HavriliakNegami(Relaxation):
|
||||
def calculation(self):
|
||||
"""Calculates the Havriliak-Negami function for
|
||||
the given parameters."""
|
||||
return self.e_inf + self.de / (1 + (1j * 2 * np.pi * self.freq * self.tau_0) ** self.alpha) ** self.beta
|
||||
return (
|
||||
self.e_inf
|
||||
+ self.de / (1 + (1j * 2 * np.pi * self.freq * self.tau_0) ** self.alpha) ** self.beta
|
||||
)
|
||||
|
||||
|
||||
class Jonscher(Relaxation):
|
||||
@@ -501,9 +567,9 @@ class Jonscher(Relaxation):
|
||||
|
||||
def calculation(self):
|
||||
"""Calculates the Q function for the given parameters"""
|
||||
return self.e_inf + (self.a_p * (2 * np.pi * self.freq / self.omega_p) ** (self.n_p - 1)) * (
|
||||
1 - 1j / np.tan(self.n_p * np.pi / 2)
|
||||
)
|
||||
return self.e_inf + (
|
||||
self.a_p * (2 * np.pi * self.freq / self.omega_p) ** (self.n_p - 1)
|
||||
) * (1 - 1j / np.tan(self.n_p * np.pi / 2))
|
||||
|
||||
|
||||
class Crim(Relaxation):
|
||||
@@ -604,7 +670,10 @@ class Crim(Relaxation):
|
||||
|
||||
def print_info(self):
|
||||
"""Print information about chosen approximation settings"""
|
||||
print(f"Approximating Complex Refractive Index Model (CRIM)" f" using {self.number_of_debye_poles} Debye poles")
|
||||
print(
|
||||
f"Approximating Complex Refractive Index Model (CRIM)"
|
||||
f" using {self.number_of_debye_poles} Debye poles"
|
||||
)
|
||||
print("CRIM parameters: ")
|
||||
for i in range(len(self.volumetric_fractions)):
|
||||
print(f"Material {i + 1}.:")
|
||||
@@ -774,17 +843,63 @@ if __name__ == "__main__":
|
||||
setup.run()
|
||||
# Testing setup
|
||||
setup = Rawdata(
|
||||
"examples/Test.txt", 0.1, 1, 0.1, "M1", number_of_debye_poles=3, plot=True, optimizer_options={"seed": 111}
|
||||
"examples/Test.txt",
|
||||
0.1,
|
||||
1,
|
||||
0.1,
|
||||
"M1",
|
||||
number_of_debye_poles=3,
|
||||
plot=True,
|
||||
optimizer_options={"seed": 111},
|
||||
)
|
||||
setup.run()
|
||||
np.random.seed(111)
|
||||
setup = HavriliakNegami(1e12, 1e-3, 0.5, 1, 10, 5, 1e-6, 0.1, 1, 0, "M2", number_of_debye_poles=6, plot=True)
|
||||
setup = HavriliakNegami(
|
||||
1e12,
|
||||
1e-3,
|
||||
0.5,
|
||||
1,
|
||||
10,
|
||||
5,
|
||||
1e-6,
|
||||
0.1,
|
||||
1,
|
||||
0,
|
||||
"M2",
|
||||
number_of_debye_poles=6,
|
||||
plot=True,
|
||||
)
|
||||
setup.run()
|
||||
setup = Jonscher(1e6, 1e-5, 50, 1, 1e5, 0.7, 0.1, 1, 0.1, "M3", number_of_debye_poles=4, plot=True)
|
||||
setup = Jonscher(
|
||||
1e6,
|
||||
1e-5,
|
||||
50,
|
||||
1,
|
||||
1e5,
|
||||
0.7,
|
||||
0.1,
|
||||
1,
|
||||
0.1,
|
||||
"M3",
|
||||
number_of_debye_poles=4,
|
||||
plot=True,
|
||||
)
|
||||
setup.run()
|
||||
f = np.array([0.5, 0.5])
|
||||
material1 = [3, 25, 1e6]
|
||||
material2 = [3, 0, 1e3]
|
||||
materials = np.array([material1, material2])
|
||||
setup = Crim(1 * 1e-1, 1e-9, 0.5, f, materials, 0.1, 1, 0, "M4", number_of_debye_poles=2, plot=True)
|
||||
setup = Crim(
|
||||
1 * 1e-1,
|
||||
1e-9,
|
||||
0.5,
|
||||
f,
|
||||
materials,
|
||||
0.1,
|
||||
1,
|
||||
0,
|
||||
"M4",
|
||||
number_of_debye_poles=2,
|
||||
plot=True,
|
||||
)
|
||||
setup.run()
|
||||
|
||||
@@ -121,7 +121,16 @@ class PSO_DLS(Optimizer):
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self, swarmsize=40, maxiter=50, omega=0.9, phip=0.9, phig=0.9, minstep=1e-8, minfun=1e-8, pflag=False, seed=None
|
||||
self,
|
||||
swarmsize=40,
|
||||
maxiter=50,
|
||||
omega=0.9,
|
||||
phip=0.9,
|
||||
phig=0.9,
|
||||
minstep=1e-8,
|
||||
minfun=1e-8,
|
||||
pflag=False,
|
||||
seed=None,
|
||||
):
|
||||
super(PSO_DLS, self).__init__(maxiter, seed)
|
||||
self.swarmsize = swarmsize
|
||||
@@ -226,10 +235,16 @@ class PSO_DLS(Optimizer):
|
||||
tmp = x[i, :].copy()
|
||||
stepsize = np.sqrt(np.sum((g - tmp) ** 2))
|
||||
if np.abs(fg - fx) <= self.minfun:
|
||||
print(f"Stopping search: Swarm best objective " f"change less than {self.minfun}")
|
||||
print(
|
||||
f"Stopping search: Swarm best objective "
|
||||
f"change less than {self.minfun}"
|
||||
)
|
||||
return tmp, fx
|
||||
elif stepsize <= self.minstep:
|
||||
print(f"Stopping search: Swarm best position " f"change less than {self.minstep}")
|
||||
print(
|
||||
f"Stopping search: Swarm best position "
|
||||
f"change less than {self.minstep}"
|
||||
)
|
||||
return tmp, fx
|
||||
else:
|
||||
g = tmp.copy()
|
||||
@@ -471,7 +486,10 @@ def DLS(logt, rl, im, freq):
|
||||
# Solving the overdetermined system y=Ax
|
||||
x = np.abs(np.linalg.lstsq(d.imag, im, rcond=None)[0])
|
||||
# x - absolute damped least-squares solution
|
||||
rp, ip = np.matmul(d.real, x[np.newaxis].T).T[0], np.matmul(d.imag, x[np.newaxis].T).T[0]
|
||||
rp, ip = (
|
||||
np.matmul(d.real, x[np.newaxis].T).T[0],
|
||||
np.matmul(d.imag, x[np.newaxis].T).T[0],
|
||||
)
|
||||
cost_i = np.sum(np.abs(ip - im)) / len(im)
|
||||
ee = np.mean(rl - rp)
|
||||
ee = max(ee, 1)
|
||||
|
||||
@@ -65,6 +65,9 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
|
||||
patchheight = 0.016
|
||||
tx = x + 0.114, y + 0.052, z + skidthickness
|
||||
else:
|
||||
logger.exception(
|
||||
"This antenna module can only be used with a spatial discretisation of 1mm or 2mm"
|
||||
)
|
||||
logger.exception(
|
||||
"This antenna module can only be used with a spatial discretisation of 1mm or 2mm"
|
||||
)
|
||||
@@ -159,13 +162,21 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
|
||||
# Metallic enclosure
|
||||
b3 = gprMax.Box(
|
||||
p1=(x + 0.025, y + casethickness, z + skidthickness),
|
||||
p2=(x + casesize[0] - 0.025, y + casesize[1] - casethickness, z + skidthickness + 0.027),
|
||||
p2=(
|
||||
x + casesize[0] - 0.025,
|
||||
y + casesize[1] - casethickness,
|
||||
z + skidthickness + 0.027,
|
||||
),
|
||||
material_id="pec",
|
||||
)
|
||||
|
||||
# Absorber material (absorber1) and foam (absorber2) around edge of absorber
|
||||
b4 = gprMax.Box(
|
||||
p1=(x + 0.025 + shieldthickness, y + casethickness + shieldthickness, z + skidthickness),
|
||||
p1=(
|
||||
x + 0.025 + shieldthickness,
|
||||
y + casethickness + shieldthickness,
|
||||
z + skidthickness,
|
||||
),
|
||||
p2=(
|
||||
x + 0.025 + shieldthickness + 0.057,
|
||||
y + casesize[1] - casethickness - shieldthickness,
|
||||
@@ -382,7 +393,11 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
|
||||
# Geometry views
|
||||
gv1 = gprMax.GeometryView(
|
||||
p1=(x - dx, y - dy, z - dz),
|
||||
p2=(x + casesize[0] + dx, y + casesize[1] + dy, z + skidthickness + casesize[2] + dz),
|
||||
p2=(
|
||||
x + casesize[0] + dx,
|
||||
y + casesize[1] + dy,
|
||||
z + skidthickness + casesize[2] + dz,
|
||||
),
|
||||
dl=(dx, dy, dz),
|
||||
filename="antenna_like_GSSI_1500",
|
||||
output_type="n",
|
||||
@@ -411,7 +426,7 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs):
|
||||
elif optstate == "GiannakisPaper":
|
||||
# Optimised custom pulse
|
||||
exc1 = gprMax.ExcitationFile(
|
||||
filepath=Path(__file__).with_name("GSSI_1500MHz_pulse.txt"),
|
||||
filepath="toolboxes/GPRAntennaModels/GSSI_1500MHz_pulse.txt",
|
||||
kind="linear",
|
||||
fill_value="extrapolate",
|
||||
)
|
||||
@@ -697,7 +712,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
|
||||
averaging="n",
|
||||
)
|
||||
b14 = gprMax.Box(
|
||||
p1=(x + 0.01 + 0.005 + 0.179, y + casethickness + 0.005 + 0.021, z + skidthickness),
|
||||
p1=(
|
||||
x + 0.01 + 0.005 + 0.179,
|
||||
y + casethickness + 0.005 + 0.021,
|
||||
z + skidthickness,
|
||||
),
|
||||
p2=(
|
||||
x + 0.01 + 0.005 + 0.195 + bowtiebase,
|
||||
y + casethickness + 0.006 + 0.202 + patchheight,
|
||||
@@ -713,7 +732,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
|
||||
# "left" side
|
||||
# extension plates
|
||||
p1 = gprMax.Plate(
|
||||
p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
|
||||
p1=(
|
||||
x + 0.01 + 0.005 + 0.025,
|
||||
y + casethickness + 0.005 + 0.021,
|
||||
z + skidthickness - 0.002,
|
||||
),
|
||||
p2=(
|
||||
x + 0.01 + 0.005 + 0.025 + bowtiebase,
|
||||
y + casethickness + 0.005 + 0.021 + patchheight,
|
||||
@@ -722,7 +745,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
|
||||
material_id="pec",
|
||||
)
|
||||
p2 = gprMax.Plate(
|
||||
p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
|
||||
p1=(
|
||||
x + 0.01 + 0.005 + 0.025,
|
||||
y + casethickness + 0.005 + 0.203,
|
||||
z + skidthickness - 0.002,
|
||||
),
|
||||
p2=(
|
||||
x + 0.01 + 0.005 + 0.025 + bowtiebase,
|
||||
y + casethickness + 0.005 + 0.203 + patchheight,
|
||||
@@ -732,7 +759,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
|
||||
)
|
||||
# triangles
|
||||
t1 = gprMax.Triangle(
|
||||
p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002),
|
||||
p1=(
|
||||
x + 0.01 + 0.005 + 0.025,
|
||||
y + casethickness + 0.005 + 0.081,
|
||||
z + skidthickness - 0.002,
|
||||
),
|
||||
p2=(
|
||||
x + 0.01 + 0.005 + 0.025 + bowtiebase,
|
||||
y + casethickness + 0.005 + 0.081,
|
||||
@@ -747,7 +778,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
|
||||
material_id="pec",
|
||||
)
|
||||
t2 = gprMax.Triangle(
|
||||
p1=(x + 0.01 + 0.005 + 0.025, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
|
||||
p1=(
|
||||
x + 0.01 + 0.005 + 0.025,
|
||||
y + casethickness + 0.005 + 0.203,
|
||||
z + skidthickness - 0.002,
|
||||
),
|
||||
p2=(
|
||||
x + 0.01 + 0.005 + 0.025 + bowtiebase,
|
||||
y + casethickness + 0.005 + 0.203,
|
||||
@@ -763,7 +798,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
|
||||
)
|
||||
# "right" side
|
||||
p3 = gprMax.Plate(
|
||||
p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002),
|
||||
p1=(
|
||||
x + 0.01 + 0.005 + 0.187,
|
||||
y + casethickness + 0.005 + 0.021,
|
||||
z + skidthickness - 0.002,
|
||||
),
|
||||
p2=(
|
||||
x + 0.01 + 0.005 + 0.187 + bowtiebase,
|
||||
y + casethickness + 0.005 + 0.021 + patchheight,
|
||||
@@ -772,7 +811,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
|
||||
material_id="pec",
|
||||
)
|
||||
p4 = gprMax.Plate(
|
||||
p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
|
||||
p1=(
|
||||
x + 0.01 + 0.005 + 0.187,
|
||||
y + casethickness + 0.005 + 0.203,
|
||||
z + skidthickness - 0.002,
|
||||
),
|
||||
p2=(
|
||||
x + 0.01 + 0.005 + 0.187 + bowtiebase,
|
||||
y + casethickness + 0.005 + 0.203 + patchheight,
|
||||
@@ -782,7 +825,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
|
||||
)
|
||||
# triangles
|
||||
t3 = gprMax.Triangle(
|
||||
p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.081, z + skidthickness - 0.002),
|
||||
p1=(
|
||||
x + 0.01 + 0.005 + 0.187,
|
||||
y + casethickness + 0.005 + 0.081,
|
||||
z + skidthickness - 0.002,
|
||||
),
|
||||
p2=(
|
||||
x + 0.01 + 0.005 + 0.187 + bowtiebase,
|
||||
y + casethickness + 0.005 + 0.081,
|
||||
@@ -797,7 +844,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
|
||||
material_id="pec",
|
||||
)
|
||||
t4 = gprMax.Triangle(
|
||||
p1=(x + 0.01 + 0.005 + 0.187, y + casethickness + 0.005 + 0.203, z + skidthickness - 0.002),
|
||||
p1=(
|
||||
x + 0.01 + 0.005 + 0.187,
|
||||
y + casethickness + 0.005 + 0.203,
|
||||
z + skidthickness - 0.002,
|
||||
),
|
||||
p2=(
|
||||
x + 0.01 + 0.005 + 0.187 + bowtiebase,
|
||||
y + casethickness + 0.005 + 0.203,
|
||||
@@ -876,7 +927,7 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
|
||||
else:
|
||||
# Optimised custom pulse
|
||||
exc1 = gprMax.ExcitationFile(
|
||||
filepath=Path(__file__).with_name("GSSI_400MHz_pulse.txt"),
|
||||
filepath="toolboxes/GPRAntennaModels/GSSI_400MHz_pulse.txt",
|
||||
kind="linear",
|
||||
fill_value="extrapolate",
|
||||
)
|
||||
@@ -906,7 +957,11 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs):
|
||||
# Geometry views
|
||||
gv1 = gprMax.GeometryView(
|
||||
p1=(x - dx, y - dy, z - dz),
|
||||
p2=(x + casesize[0] + dx, y + casesize[1] + dy, z + skidthickness + casesize[2] + dz),
|
||||
p2=(
|
||||
x + casesize[0] + dx,
|
||||
y + casesize[1] + dy,
|
||||
z + skidthickness + casesize[2] + dz,
|
||||
),
|
||||
dl=(dx, dy, dz),
|
||||
filename="antenna_like_GSSI_400",
|
||||
output_type="n",
|
||||
|
||||
某些文件未显示，因为此 diff 中更改的文件太多 显示更多
在新工单中引用
屏蔽一个用户