diff --git a/CREDITS b/CREDITS index 506c8bf2..bdf69641 100644 --- a/CREDITS +++ b/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 diff --git a/docs/.gitignore b/docs/.gitignore new file mode 100644 index 00000000..e775aef0 --- /dev/null +++ b/docs/.gitignore @@ -0,0 +1,8 @@ +_*/ +doctrees/ +dirhtml/ + +*.html +.buildinfo +objects.inv +searchindex.js diff --git a/docs/source/input_api.rst b/docs/source/input_api.rst index 94f6790e..a64ba49a 100644 --- a/docs/source/input_api.rst +++ b/docs/source/input_api.rst @@ -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 diff --git a/docs/source/input_hash_cmds.rst b/docs/source/input_hash_cmds.rst index 47002aba..a8ed9699 100644 --- a/docs/source/input_hash_cmds.rst +++ b/docs/source/input_hash_cmds.rst @@ -622,7 +622,7 @@ For example, to create an orthogonal parallelepiped with fractal distributed pro .. note:: - * Currently (2024) we are not aware of a formulation of Perfectly Matched Layer (PML) absorbing boundary that can specifically handle distributions of material properties (such as those created by fractals) throughout the thickness of the PML, i.e. this is a required area of research. Our PML formulations can work to an extent depending on your modelling scenario and requirements. You may need to increase the thickness of the PML and/or consider tuning the parameters of the PML (:ref:`pml-tuning`) to improve performance for your specific model. + * Currently (2024) we are not aware of a formulation of Perfectly Matched Layer (PML) absorbing boundary that can specifically handle distributions of material properties (such as those created by fractals) throughout the thickness of the PML, i.e. this is a required area of research. Our PML formulations can work to an extent depending on your modelling scenario and requirements. You may need to increase the thickness of the PML and/or consider tuning the parameters of the PML (:ref:`pml-tuning`) to improve performance for your specific model. #add_surface_roughness: ----------------------- @@ -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 ` 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 `_ (GSoC) project and `more details can be found in the original pull request `_. + #rx: ---- diff --git a/docs/source/references.rst b/docs/source/references.rst index 954c1b04..b416a69b 100644 --- a/docs/source/references.rst +++ b/docs/source/references.rst @@ -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) diff --git a/gprMax/__init__.py b/gprMax/__init__.py index 0dba095a..1b35175e 100644 --- a/gprMax/__init__.py +++ b/gprMax/__init__.py @@ -32,6 +32,7 @@ from .user_objects.cmds_multiuse import ( AddDebyeDispersion, AddDrudeDispersion, AddLorentzDispersion, + DiscretePlaneWave, ExcitationFile, HertzianDipole, MagneticDipole, diff --git a/gprMax/config.pxd b/gprMax/config.pxd index cdd76ac7..a30ce6a0 100644 --- a/gprMax/config.pxd +++ b/gprMax/config.pxd @@ -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. diff --git a/gprMax/config.py b/gprMax/config.py index 06b383ac..67f9d328 100644 --- a/gprMax/config.py +++ b/gprMax/config.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. @@ -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. diff --git a/gprMax/contexts.py b/gprMax/contexts.py index dc956da9..57e980d2 100644 --- a/gprMax/contexts.py +++ b/gprMax/contexts.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. diff --git a/gprMax/cuda_opencl/knl_fields_updates.py b/gprMax/cuda_opencl/knl_fields_updates.py index f6bd614d..3784595b 100644 --- a/gprMax/cuda_opencl/knl_fields_updates.py +++ b/gprMax/cuda_opencl/knl_fields_updates.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. @@ -238,7 +238,7 @@ update_electric_dispersive_A = { int x_T = (i % ($NX_T * $NY_T * $NZ_T)) / ($NY_T * $NZ_T); int y_T = ((i % ($NX_T * $NY_T * $NZ_T)) % ($NY_T * $NZ_T)) / $NZ_T; int z_T = ((i % ($NX_T * $NY_T * $NZ_T)) % ($NY_T * $NZ_T)) % $NZ_T; - + // Ex component if ((NY != 1 || NZ != 1) && x >= 0 && x < NX && y > 0 && y < NY && z > 0 && z < NZ) { int materialEx = ID[IDX4D_ID(0,x_ID,y_ID,z_ID)]; diff --git a/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py b/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py index 7382504c..de17c785 100644 --- a/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py +++ b/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.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. diff --git a/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py b/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py index e56788ee..59744af5 100644 --- a/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py +++ b/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.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. diff --git a/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py b/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py index c16b48ff..e8d35feb 100644 --- a/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py +++ b/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.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. diff --git a/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py b/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py index 384738d4..ea79c347 100644 --- a/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py +++ b/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.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. diff --git a/gprMax/cuda_opencl/knl_snapshots.py b/gprMax/cuda_opencl/knl_snapshots.py index 07addd30..d7a799d8 100644 --- a/gprMax/cuda_opencl/knl_snapshots.py +++ b/gprMax/cuda_opencl/knl_snapshots.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. diff --git a/gprMax/cuda_opencl/knl_source_updates.py b/gprMax/cuda_opencl/knl_source_updates.py index 8ee43ecb..37932a8c 100644 --- a/gprMax/cuda_opencl/knl_source_updates.py +++ b/gprMax/cuda_opencl/knl_source_updates.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. diff --git a/gprMax/cuda_opencl/knl_store_outputs.py b/gprMax/cuda_opencl/knl_store_outputs.py index 85371903..54e844db 100644 --- a/gprMax/cuda_opencl/knl_store_outputs.py +++ b/gprMax/cuda_opencl/knl_store_outputs.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. diff --git a/gprMax/cython/fields_updates_dispersive_template.jinja b/gprMax/cython/fields_updates_dispersive_template.jinja index eac7601a..dff865a1 100644 --- a/gprMax/cython/fields_updates_dispersive_template.jinja +++ b/gprMax/cython/fields_updates_dispersive_template.jinja @@ -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. diff --git a/gprMax/cython/fields_updates_hsg.pyx b/gprMax/cython/fields_updates_hsg.pyx index 878999ee..d5c72c04 100644 --- a/gprMax/cython/fields_updates_hsg.pyx +++ b/gprMax/cython/fields_updates_hsg.pyx @@ -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. diff --git a/gprMax/cython/fields_updates_normal.pyx b/gprMax/cython/fields_updates_normal.pyx index 6eb8fae9..035bd0b2 100644 --- a/gprMax/cython/fields_updates_normal.pyx +++ b/gprMax/cython/fields_updates_normal.pyx @@ -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. diff --git a/gprMax/cython/fractals_generate.pyx b/gprMax/cython/fractals_generate.pyx index a3f101b5..2339b6bf 100644 --- a/gprMax/cython/fractals_generate.pyx +++ b/gprMax/cython/fractals_generate.pyx @@ -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. diff --git a/gprMax/cython/geometry_outputs.pyx b/gprMax/cython/geometry_outputs.pyx index df074aee..15006308 100644 --- a/gprMax/cython/geometry_outputs.pyx +++ b/gprMax/cython/geometry_outputs.pyx @@ -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. diff --git a/gprMax/cython/geometry_primitives.pyx b/gprMax/cython/geometry_primitives.pyx index 89d3a92d..4d1900ca 100644 --- a/gprMax/cython/geometry_primitives.pyx +++ b/gprMax/cython/geometry_primitives.pyx @@ -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: diff --git a/gprMax/cython/plane_wave.pyx b/gprMax/cython/plane_wave.pyx new file mode 100644 index 00000000..408f59b6 --- /dev/null +++ b/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 . + +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]/integers[2, 1]-required_ratio)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] 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) diff --git a/gprMax/cython/pml_build.pyx b/gprMax/cython/pml_build.pyx index 9299207c..fb369db4 100644 --- a/gprMax/cython/pml_build.pyx +++ b/gprMax/cython/pml_build.pyx @@ -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. diff --git a/gprMax/cython/pml_updates_electric_HORIPML.pyx b/gprMax/cython/pml_updates_electric_HORIPML.pyx index f73530c4..19571c65 100644 --- a/gprMax/cython/pml_updates_electric_HORIPML.pyx +++ b/gprMax/cython/pml_updates_electric_HORIPML.pyx @@ -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. diff --git a/gprMax/cython/pml_updates_electric_MRIPML.pyx b/gprMax/cython/pml_updates_electric_MRIPML.pyx index 44d497cf..e79de36d 100644 --- a/gprMax/cython/pml_updates_electric_MRIPML.pyx +++ b/gprMax/cython/pml_updates_electric_MRIPML.pyx @@ -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. diff --git a/gprMax/cython/pml_updates_magnetic_HORIPML.pyx b/gprMax/cython/pml_updates_magnetic_HORIPML.pyx index 1b0ea55a..7843d346 100644 --- a/gprMax/cython/pml_updates_magnetic_HORIPML.pyx +++ b/gprMax/cython/pml_updates_magnetic_HORIPML.pyx @@ -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. diff --git a/gprMax/cython/pml_updates_magnetic_MRIPML.pyx b/gprMax/cython/pml_updates_magnetic_MRIPML.pyx index fd4ac73b..02934497 100644 --- a/gprMax/cython/pml_updates_magnetic_MRIPML.pyx +++ b/gprMax/cython/pml_updates_magnetic_MRIPML.pyx @@ -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. diff --git a/gprMax/cython/snapshots.pyx b/gprMax/cython/snapshots.pyx index ec91e3b8..c5b26833 100644 --- a/gprMax/cython/snapshots.pyx +++ b/gprMax/cython/snapshots.pyx @@ -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. diff --git a/gprMax/cython/yee_cell_build.pyx b/gprMax/cython/yee_cell_build.pyx index 1412561f..4814fcd5 100644 --- a/gprMax/cython/yee_cell_build.pyx +++ b/gprMax/cython/yee_cell_build.pyx @@ -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. diff --git a/gprMax/cython/yee_cell_setget_rigid.pxd b/gprMax/cython/yee_cell_setget_rigid.pxd index 001124ce..4461ade2 100644 --- a/gprMax/cython/yee_cell_setget_rigid.pxd +++ b/gprMax/cython/yee_cell_setget_rigid.pxd @@ -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. diff --git a/gprMax/cython/yee_cell_setget_rigid.pyx b/gprMax/cython/yee_cell_setget_rigid.pyx index 99fb470d..ff36e9db 100644 --- a/gprMax/cython/yee_cell_setget_rigid.pyx +++ b/gprMax/cython/yee_cell_setget_rigid.pyx @@ -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. diff --git a/gprMax/fields_outputs.py b/gprMax/fields_outputs.py index ce4b9197..3365b7ba 100644 --- a/gprMax/fields_outputs.py +++ b/gprMax/fields_outputs.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. @@ -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 diff --git a/gprMax/fractals.py b/gprMax/fractals.py index d2670a1f..bfe9e5be 100644 --- a/gprMax/fractals.py +++ b/gprMax/fractals.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. @@ -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 diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py index 9f971568..f48d5110 100644 --- a/gprMax/geometry_outputs.py +++ b/gprMax/geometry_outputs.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. diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index e214ac5d..34664071 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.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. diff --git a/gprMax/grid/fdtd_grid.py b/gprMax/grid/fdtd_grid.py index 265bee08..30306a67 100644 --- a/gprMax/grid/fdtd_grid.py +++ b/gprMax/grid/fdtd_grid.py @@ -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] diff --git a/gprMax/hash_cmds_file.py b/gprMax/hash_cmds_file.py index 283cb50c..b177c175 100644 --- a/gprMax/hash_cmds_file.py +++ b/gprMax/hash_cmds_file.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. @@ -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", diff --git a/gprMax/hash_cmds_geometry.py b/gprMax/hash_cmds_geometry.py index 1cb068d3..81e4476b 100644 --- a/gprMax/hash_cmds_geometry.py +++ b/gprMax/hash_cmds_geometry.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. @@ -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: diff --git a/gprMax/hash_cmds_multiuse.py b/gprMax/hash_cmds_multiuse.py index 9e6d0b55..97040747 100644 --- a/gprMax/hash_cmds_multiuse.py +++ b/gprMax/hash_cmds_multiuse.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. @@ -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) diff --git a/gprMax/hash_cmds_singleuse.py b/gprMax/hash_cmds_singleuse.py index 08b57b97..fd6303fd 100644 --- a/gprMax/hash_cmds_singleuse.py +++ b/gprMax/hash_cmds_singleuse.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. diff --git a/gprMax/materials.py b/gprMax/materials.py index 05c35af0..71ed23e9 100644 --- a/gprMax/materials.py +++ b/gprMax/materials.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. diff --git a/gprMax/pml.py b/gprMax/pml.py index 8b94fb08..6bfe9dbf 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.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. @@ -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 diff --git a/gprMax/receivers.py b/gprMax/receivers.py index f27dae23..248966f4 100644 --- a/gprMax/receivers.py +++ b/gprMax/receivers.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. diff --git a/gprMax/scene.py b/gprMax/scene.py index df7ae687..5af3d29a 100644 --- a/gprMax/scene.py +++ b/gprMax/scene.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. diff --git a/gprMax/snapshots.py b/gprMax/snapshots.py index c695c7fe..60b3722d 100644 --- a/gprMax/snapshots.py +++ b/gprMax/snapshots.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. diff --git a/gprMax/solvers.py b/gprMax/solvers.py index eadfdbe8..a3d58544 100644 --- a/gprMax/solvers.py +++ b/gprMax/solvers.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. @@ -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): diff --git a/gprMax/sources.py b/gprMax/sources.py old mode 100644 new mode 100755 index 1ed6782d..e5f7d946 --- a/gprMax/sources.py +++ b/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 + ) diff --git a/gprMax/subgrids/grid.py b/gprMax/subgrids/grid.py index aa33d958..66e7f02f 100644 --- a/gprMax/subgrids/grid.py +++ b/gprMax/subgrids/grid.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. diff --git a/gprMax/subgrids/precursor_nodes.py b/gprMax/subgrids/precursor_nodes.py index cd421fef..b162c457 100644 --- a/gprMax/subgrids/precursor_nodes.py +++ b/gprMax/subgrids/precursor_nodes.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. @@ -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: diff --git a/gprMax/subgrids/subgrid_hsg.py b/gprMax/subgrids/subgrid_hsg.py index 62e50fc7..194f48cb 100644 --- a/gprMax/subgrids/subgrid_hsg.py +++ b/gprMax/subgrids/subgrid_hsg.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. @@ -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)" diff --git a/gprMax/subgrids/updates.py b/gprMax/subgrids/updates.py index 266d87c0..8a0a04d7 100644 --- a/gprMax/subgrids/updates.py +++ b/gprMax/subgrids/updates.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. diff --git a/gprMax/subgrids/user_objects.py b/gprMax/subgrids/user_objects.py index b6d1e02b..c306b9fe 100644 --- a/gprMax/subgrids/user_objects.py +++ b/gprMax/subgrids/user_objects.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. diff --git a/gprMax/taskfarm.py b/gprMax/taskfarm.py index 84f779d5..335a3ae1 100644 --- a/gprMax/taskfarm.py +++ b/gprMax/taskfarm.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: Tobias Schruff # # This file is part of gprMax. diff --git a/gprMax/updates/cpu_updates.py b/gprMax/updates/cpu_updates.py index 8feecc88..b4a9ef32 100644 --- a/gprMax/updates/cpu_updates.py +++ b/gprMax/updates/cpu_updates.py @@ -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. diff --git a/gprMax/user_inputs.py b/gprMax/user_inputs.py index 2b7fd94c..3ac9a5ad 100644 --- a/gprMax/user_inputs.py +++ b/gprMax/user_inputs.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. diff --git a/gprMax/user_objects/cmds_geometry/add_grass.py b/gprMax/user_objects/cmds_geometry/add_grass.py index 87a0958c..57345e4c 100644 --- a/gprMax/user_objects/cmds_geometry/add_grass.py +++ b/gprMax/user_objects/cmds_geometry/add_grass.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. @@ -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( diff --git a/gprMax/user_objects/cmds_geometry/add_surface_roughness.py b/gprMax/user_objects/cmds_geometry/add_surface_roughness.py index fc359f13..5fe591cb 100644 --- a/gprMax/user_objects/cmds_geometry/add_surface_roughness.py +++ b/gprMax/user_objects/cmds_geometry/add_surface_roughness.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. @@ -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: diff --git a/gprMax/user_objects/cmds_geometry/add_surface_water.py b/gprMax/user_objects/cmds_geometry/add_surface_water.py index 1d7b2b23..72e5b3d2 100644 --- a/gprMax/user_objects/cmds_geometry/add_surface_water.py +++ b/gprMax/user_objects/cmds_geometry/add_surface_water.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. diff --git a/gprMax/user_objects/cmds_geometry/box.py b/gprMax/user_objects/cmds_geometry/box.py index 4c715049..ef33a5c7 100644 --- a/gprMax/user_objects/cmds_geometry/box.py +++ b/gprMax/user_objects/cmds_geometry/box.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. @@ -147,7 +147,6 @@ class Box(RotatableMixin, GeometryUserObject): yf, zs, zf, - config.get_model_config().ompthreads, numID, numIDx, numIDy, diff --git a/gprMax/user_objects/cmds_geometry/build_templates.py b/gprMax/user_objects/cmds_geometry/build_templates.py index e3298917..84d17f66 100644 --- a/gprMax/user_objects/cmds_geometry/build_templates.py +++ b/gprMax/user_objects/cmds_geometry/build_templates.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. diff --git a/gprMax/user_objects/cmds_geometry/cmds_geometry.py b/gprMax/user_objects/cmds_geometry/cmds_geometry.py index 13af63dd..df624c4a 100644 --- a/gprMax/user_objects/cmds_geometry/cmds_geometry.py +++ b/gprMax/user_objects/cmds_geometry/cmds_geometry.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. diff --git a/gprMax/user_objects/cmds_geometry/cone.py b/gprMax/user_objects/cmds_geometry/cone.py index fffa1729..60b269ef 100644 --- a/gprMax/user_objects/cmds_geometry/cone.py +++ b/gprMax/user_objects/cmds_geometry/cone.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. diff --git a/gprMax/user_objects/cmds_geometry/cylinder.py b/gprMax/user_objects/cmds_geometry/cylinder.py index 7bcb6c15..0904f039 100644 --- a/gprMax/user_objects/cmds_geometry/cylinder.py +++ b/gprMax/user_objects/cmds_geometry/cylinder.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. diff --git a/gprMax/user_objects/cmds_geometry/cylindrical_sector.py b/gprMax/user_objects/cmds_geometry/cylindrical_sector.py index 97da02ca..af2cb4b4 100644 --- a/gprMax/user_objects/cmds_geometry/cylindrical_sector.py +++ b/gprMax/user_objects/cmds_geometry/cylindrical_sector.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. diff --git a/gprMax/user_objects/cmds_geometry/edge.py b/gprMax/user_objects/cmds_geometry/edge.py index 9922b774..69edd75b 100644 --- a/gprMax/user_objects/cmds_geometry/edge.py +++ b/gprMax/user_objects/cmds_geometry/edge.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. diff --git a/gprMax/user_objects/cmds_geometry/ellipsoid.py b/gprMax/user_objects/cmds_geometry/ellipsoid.py index 464c351f..f38911f4 100644 --- a/gprMax/user_objects/cmds_geometry/ellipsoid.py +++ b/gprMax/user_objects/cmds_geometry/ellipsoid.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. diff --git a/gprMax/user_objects/cmds_geometry/fractal_box.py b/gprMax/user_objects/cmds_geometry/fractal_box.py index 48a22dbf..f38ad768 100644 --- a/gprMax/user_objects/cmds_geometry/fractal_box.py +++ b/gprMax/user_objects/cmds_geometry/fractal_box.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. @@ -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, diff --git a/gprMax/user_objects/cmds_geometry/geometry_objects_read.py b/gprMax/user_objects/cmds_geometry/geometry_objects_read.py index 789a569e..53b41105 100644 --- a/gprMax/user_objects/cmds_geometry/geometry_objects_read.py +++ b/gprMax/user_objects/cmds_geometry/geometry_objects_read.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. @@ -151,7 +151,6 @@ class GeometryObjectsRead(GeometryUserObject): xs, ys, zs, - config.get_model_config().ompthreads, numexistmaterials, averaging, data, diff --git a/gprMax/user_objects/cmds_geometry/plate.py b/gprMax/user_objects/cmds_geometry/plate.py index 6da03d1a..a2c90195 100644 --- a/gprMax/user_objects/cmds_geometry/plate.py +++ b/gprMax/user_objects/cmds_geometry/plate.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. diff --git a/gprMax/user_objects/cmds_geometry/sphere.py b/gprMax/user_objects/cmds_geometry/sphere.py index ef2db9cd..33364066 100644 --- a/gprMax/user_objects/cmds_geometry/sphere.py +++ b/gprMax/user_objects/cmds_geometry/sphere.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. diff --git a/gprMax/user_objects/cmds_geometry/triangle.py b/gprMax/user_objects/cmds_geometry/triangle.py index 27d55edc..54fc3ae5 100644 --- a/gprMax/user_objects/cmds_geometry/triangle.py +++ b/gprMax/user_objects/cmds_geometry/triangle.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. diff --git a/gprMax/user_objects/cmds_multiuse.py b/gprMax/user_objects/cmds_multiuse.py index ae019f9f..49891726 100644 --- a/gprMax/user_objects/cmds_multiuse.py +++ b/gprMax/user_objects/cmds_multiuse.py @@ -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. diff --git a/gprMax/utilities/host_info.py b/gprMax/utilities/host_info.py index 482c099d..e7d2aa4c 100644 --- a/gprMax/utilities/host_info.py +++ b/gprMax/utilities/host_info.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. @@ -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 diff --git a/gprMax/utilities/logging.py b/gprMax/utilities/logging.py index 3760698f..aa1a6640 100644 --- a/gprMax/utilities/logging.py +++ b/gprMax/utilities/logging.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. diff --git a/gprMax/utilities/utilities.py b/gprMax/utilities/utilities.py index 74ac6dc4..f28e53da 100644 --- a/gprMax/utilities/utilities.py +++ b/gprMax/utilities/utilities.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. @@ -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 __ __ __ _ _ __ _ __| \/ | __ ___ __ / _` | '_ \| '__| |\/| |/ _` \ \/ / | (_| | |_) | | | | | | (_| |> < diff --git a/gprMax/waveforms.py b/gprMax/waveforms.py index 6c2de0ba..19f02f7a 100644 --- a/gprMax/waveforms.py +++ b/gprMax/waveforms.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. @@ -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": diff --git a/gprMax.toml b/pyproject.toml similarity index 65% rename from gprMax.toml rename to pyproject.toml index ae219ea9..ac7fe4f4 100644 --- a/gprMax.toml +++ b/pyproject.toml @@ -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" diff --git a/reframe_tests/regression_checks/TestDiscretePlaneWaveSource_a0061267/discrete_plane_wave_fs.h5 b/reframe_tests/regression_checks/TestDiscretePlaneWaveSource_a0061267/discrete_plane_wave_fs.h5 new file mode 100644 index 00000000..b6bfa97c Binary files /dev/null and b/reframe_tests/regression_checks/TestDiscretePlaneWaveSource_a0061267/discrete_plane_wave_fs.h5 differ diff --git a/reframe_tests/regression_checks/TestEdgeGeometry_c4eb9ee8/antenna_wire_dipole_fs.h5 b/reframe_tests/regression_checks/TestEdgeGeometry_c4eb9ee8/antenna_wire_dipole_fs.h5 new file mode 100644 index 00000000..a1eb5455 Binary files /dev/null and b/reframe_tests/regression_checks/TestEdgeGeometry_c4eb9ee8/antenna_wire_dipole_fs.h5 differ diff --git a/reframe_tests/tests/src/source_tests/discrete_plane_wave_fs.in b/reframe_tests/tests/src/source_tests/discrete_plane_wave_fs.in new file mode 100644 index 00000000..0c8cc4db --- /dev/null +++ b/reframe_tests/tests/src/source_tests/discrete_plane_wave_fs.in @@ -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 diff --git a/reframe_tests/tests/test_sources.py b/reframe_tests/tests/test_sources.py index 9bca3b42..5cc9f216 100644 --- a/reframe_tests/tests/test_sources.py +++ b/reframe_tests/tests/test_sources.py @@ -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"} diff --git a/setup.py b/setup.py index 409a97bd..456c0683 100644 --- a/setup.py +++ b/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+", diff --git a/testing/analytical_solutions.py b/testing/analytical_solutions.py index bbbd5f94..94256f94 100644 --- a/testing/analytical_solutions.py +++ b/testing/analytical_solutions.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. @@ -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 diff --git a/testing/benchmarking/bench_simple.py b/testing/benchmarking/bench_simple.py index 7898aaa2..0a3bd6e7 100644 --- a/testing/benchmarking/bench_simple.py +++ b/testing/benchmarking/bench_simple.py @@ -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 diff --git a/testing/diff_output_files.py b/testing/diff_output_files.py index 0ee81d36..46c2f1a7 100644 --- a/testing/diff_output_files.py +++ b/testing/diff_output_files.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. @@ -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) diff --git a/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py b/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py index f958259d..6a60cbf0 100644 --- a/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py +++ b/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.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. @@ -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)) diff --git a/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate.py b/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate.py index cac3b3e3..3a352903 100755 --- a/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate.py +++ b/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate.py @@ -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") diff --git a/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate_ref.py b/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate_ref.py index 97e2d75c..80658143 100755 --- a/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate_ref.py +++ b/testing/models_pmls/pml_3D_pec_plate/pml_3D_pec_plate_ref.py @@ -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") diff --git a/testing/models_pmls/pml_basic/plot_pml_comparison.py b/testing/models_pmls/pml_basic/plot_pml_comparison.py index 767ca972..6a04d542 100644 --- a/testing/models_pmls/pml_basic/plot_pml_comparison.py +++ b/testing/models_pmls/pml_basic/plot_pml_comparison.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. @@ -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() diff --git a/testing/models_pmls/pml_basic/pml_basic.py b/testing/models_pmls/pml_basic/pml_basic.py index 0c29b4f1..a54492e8 100755 --- a/testing/models_pmls/pml_basic/pml_basic.py +++ b/testing/models_pmls/pml_basic/pml_basic.py @@ -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 diff --git a/testing/other_codes/vs_MoM_MATLAB/antenna_bowtie_fs/antenna_bowtie_fs.py b/testing/other_codes/vs_MoM_MATLAB/antenna_bowtie_fs/antenna_bowtie_fs.py index fd2cdbfc..8d1d666f 100644 --- a/testing/other_codes/vs_MoM_MATLAB/antenna_bowtie_fs/antenna_bowtie_fs.py +++ b/testing/other_codes/vs_MoM_MATLAB/antenna_bowtie_fs/antenna_bowtie_fs.py @@ -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", diff --git a/testing/test_experimental.py b/testing/test_experimental.py index 62118c5e..b8b7921b 100644 --- a/testing/test_experimental.py +++ b/testing/test_experimental.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. @@ -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") diff --git a/testing/test_models.py b/testing/test_models.py index 10241217..4afa8ae8 100644 --- a/testing/test_models.py +++ b/testing/test_models.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. @@ -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( diff --git a/toolboxes/AntennaPatterns/plot_fields.py b/toolboxes/AntennaPatterns/plot_fields.py index b43160ed..f79aeed7 100644 --- a/toolboxes/AntennaPatterns/plot_fields.py +++ b/toolboxes/AntennaPatterns/plot_fields.py @@ -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] diff --git a/toolboxes/AustinManWoman/head_only_h5.py b/toolboxes/AustinManWoman/head_only_h5.py index ff4de7fb..0ff841dc 100644 --- a/toolboxes/AustinManWoman/head_only_h5.py +++ b/toolboxes/AustinManWoman/head_only_h5.py @@ -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" diff --git a/toolboxes/DebyeFit/Debye_Fit.py b/toolboxes/DebyeFit/Debye_Fit.py index 47fa6aab..37d76133 100644 --- a/toolboxes/DebyeFit/Debye_Fit.py +++ b/toolboxes/DebyeFit/Debye_Fit.py @@ -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() diff --git a/toolboxes/DebyeFit/optimization.py b/toolboxes/DebyeFit/optimization.py index 0953bb6f..aeb7422f 100644 --- a/toolboxes/DebyeFit/optimization.py +++ b/toolboxes/DebyeFit/optimization.py @@ -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) diff --git a/toolboxes/GPRAntennaModels/GSSI.py b/toolboxes/GPRAntennaModels/GSSI.py index 582b839d..4487a224 100644 --- a/toolboxes/GPRAntennaModels/GSSI.py +++ b/toolboxes/GPRAntennaModels/GSSI.py @@ -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 --git a/toolboxes/GPRAntennaModels/MALA.py b/toolboxes/GPRAntennaModels/MALA.py index cbee5f06..585822a2 100644 --- a/toolboxes/GPRAntennaModels/MALA.py +++ b/toolboxes/GPRAntennaModels/MALA.py @@ -84,20 +84,30 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): bowtieheight = 0.024 tx = x + 0.062, y + 0.052, z + skidthickness else: - logger.exception("This antenna module can only be used with a spatial resolution of 1mm or 2mm") + logger.exception( + "This antenna module can only be used with a spatial resolution of 1mm or 2mm" + ) raise ValueError # SMD resistors - 3 on each Tx & Rx bowtie arm txres = 470 # Ohms txrescellupper = txres / 3 # Resistor over 3 cells - txsigupper = ((1 / txrescellupper) * (dy / (dx * dz))) / 2 # Divide by number of parallel edges per resistor + txsigupper = ( + (1 / txrescellupper) * (dy / (dx * dz)) + ) / 2 # Divide by number of parallel edges per resistor txrescelllower = txres / 4 # Resistor over 4 cells - txsiglower = ((1 / txrescelllower) * (dy / (dx * dz))) / 2 # Divide by number of parallel edges per resistor + txsiglower = ( + (1 / txrescelllower) * (dy / (dx * dz)) + ) / 2 # Divide by number of parallel edges per resistor rxres = 150 # Ohms rxrescellupper = rxres / 3 # Resistor over 3 cells - rxsigupper = ((1 / rxrescellupper) * (dy / (dx * dz))) / 2 # Divide by number of parallel edges per resistor + rxsigupper = ( + (1 / rxrescellupper) * (dy / (dx * dz)) + ) / 2 # Divide by number of parallel edges per resistor rxrescelllower = rxres / 4 # Resistor over 4 cells - rxsiglower = ((1 / rxrescelllower) * (dy / (dx * dz))) / 2 # Divide by number of parallel edges per resistor + rxsiglower = ( + (1 / rxrescelllower) * (dy / (dx * dz)) + ) / 2 # Divide by number of parallel edges per resistor # Material definitions absorber = gprMax.Material(er=absorberEr, se=absorbersig, mr=1, sm=0, id="absorber") @@ -108,7 +118,18 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): txresupper = gprMax.Material(er=3, se=txsigupper, mr=1, sm=0, id="txresupper") rxreslower = gprMax.Material(er=3, se=rxsiglower, mr=1, sm=0, id="rxreslower") rxresupper = gprMax.Material(er=3, se=rxsigupper, mr=1, sm=0, id="rxresupper") - scene_objects.extend((absorber, pcb, hdpe, polypropylene, txreslower, txresupper, rxreslower, rxresupper)) + scene_objects.extend( + ( + absorber, + pcb, + hdpe, + polypropylene, + txreslower, + txresupper, + rxreslower, + rxresupper, + ) + ) # Antenna geometry # Shield - metallic enclosure @@ -119,19 +140,31 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): ) b2 = gprMax.Box( p1=(x + 0.020, y + casethickness, z + skidthickness), - p2=(x + 0.100, y + casesize[1] - casethickness, z + skidthickness + casethickness), + p2=( + x + 0.100, + y + casesize[1] - casethickness, + z + skidthickness + casethickness, + ), material_id="free_space", ) b3 = gprMax.Box( p1=(x + 0.100, y + casethickness, z + skidthickness), - p2=(x + casesize[0] - casethickness, y + casesize[1] - casethickness, z + skidthickness + casethickness), + p2=( + x + casesize[0] - casethickness, + y + casesize[1] - casethickness, + z + skidthickness + casethickness, + ), material_id="free_space", ) # Absorber material b4 = gprMax.Box( p1=(x + 0.020, y + casethickness, z + skidthickness), - p2=(x + 0.100, y + casesize[1] - casethickness, z + skidthickness + casesize[2] - casethickness), + p2=( + x + 0.100, + y + casesize[1] - casethickness, + z + skidthickness + casesize[2] - casethickness, + ), material_id="absorber", ) b5 = gprMax.Box( @@ -148,7 +181,11 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): # Shield - cylindrical sections c1 = gprMax.Cylinder( p1=(x + 0.055, y + casesize[1] - 0.008, z + skidthickness), - p2=(x + 0.055, y + casesize[1] - 0.008, z + skidthickness + casesize[2] - casethickness), + p2=( + x + 0.055, + y + casesize[1] - 0.008, + z + skidthickness + casesize[2] - casethickness, + ), r=0.008, material_id="pec", ) @@ -160,7 +197,11 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): ) c3 = gprMax.Cylinder( p1=(x + 0.147, y + casesize[1] - 0.008, z + skidthickness), - p2=(x + 0.147, y + casesize[1] - 0.008, z + skidthickness + casesize[2] - casethickness), + p2=( + x + 0.147, + y + casesize[1] - 0.008, + z + skidthickness + casesize[2] - casethickness, + ), r=0.008, material_id="pec", ) @@ -172,7 +213,11 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): ) c5 = gprMax.Cylinder( p1=(x + 0.055, y + casesize[1] - 0.008, z + skidthickness), - p2=(x + 0.055, y + casesize[1] - 0.008, z + skidthickness + casesize[2] - casethickness), + p2=( + x + 0.055, + y + casesize[1] - 0.008, + z + skidthickness + casesize[2] - casethickness, + ), r=0.007, material_id="free_space", ) @@ -184,7 +229,11 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): ) c7 = gprMax.Cylinder( p1=(x + 0.147, y + casesize[1] - 0.008, z + skidthickness), - p2=(x + 0.147, y + casesize[1] - 0.008, z + skidthickness + casesize[2] - casethickness), + p2=( + x + 0.147, + y + casesize[1] - 0.008, + z + skidthickness + casesize[2] - casethickness, + ), r=0.007, material_id="free_space", ) @@ -196,7 +245,11 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): ) b6 = gprMax.Box( p1=(x + 0.054, y + casesize[1] - 0.016, z + skidthickness), - p2=(x + 0.056, y + casesize[1] - 0.014, z + skidthickness + casesize[2] - casethickness), + p2=( + x + 0.056, + y + casesize[1] - 0.014, + z + skidthickness + casesize[2] - casethickness, + ), material_id="free_space", ) b7 = gprMax.Box( @@ -206,7 +259,11 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): ) b8 = gprMax.Box( p1=(x + 0.146, y + casesize[1] - 0.016, z + skidthickness), - p2=(x + 0.148, y + casesize[1] - 0.014, z + skidthickness + casesize[2] - casethickness), + p2=( + x + 0.148, + y + casesize[1] - 0.014, + z + skidthickness + casesize[2] - casethickness, + ), material_id="free_space", ) b9 = gprMax.Box( @@ -219,18 +276,30 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): # PCB b10 = gprMax.Box( p1=(x + 0.020, y + 0.018, z + skidthickness), - p2=(x + casesize[0] - casethickness, y + casesize[1] - 0.018, z + skidthickness + pcbthickness), + p2=( + x + casesize[0] - casethickness, + y + casesize[1] - 0.018, + z + skidthickness + pcbthickness, + ), material_id="pcb", ) # Shield - Tx & Rx cavities b11 = gprMax.Box( p1=(x + 0.032, y + 0.022, z + skidthickness), - p2=(x + 0.032 + cavitysize[0], y + 0.022 + cavitysize[1], z + skidthickness + cavitysize[2]), + p2=( + x + 0.032 + cavitysize[0], + y + 0.022 + cavitysize[1], + z + skidthickness + cavitysize[2], + ), material_id="pec", ) b12 = gprMax.Box( - p1=(x + 0.032 + cavitythickness, y + 0.022 + cavitythickness, z + skidthickness), + p1=( + x + 0.032 + cavitythickness, + y + 0.022 + cavitythickness, + z + skidthickness, + ), p2=( x + 0.032 + cavitysize[0] - cavitythickness, y + 0.022 + cavitysize[1] - cavitythickness, @@ -240,11 +309,19 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): ) b13 = gprMax.Box( p1=(x + 0.108, y + 0.022, z + skidthickness), - p2=(x + 0.108 + cavitysize[0], y + 0.022 + cavitysize[1], z + skidthickness + cavitysize[2]), + p2=( + x + 0.108 + cavitysize[0], + y + 0.022 + cavitysize[1], + z + skidthickness + cavitysize[2], + ), material_id="pec", ) b14 = gprMax.Box( - p1=(x + 0.108 + cavitythickness, y + 0.022 + cavitythickness, z + skidthickness), + p1=( + x + 0.108 + cavitythickness, + y + 0.022 + cavitythickness, + z + skidthickness, + ), p2=( x + 0.108 + cavitysize[0] - cavitythickness, y + 0.022 + cavitysize[1] - cavitythickness, @@ -264,14 +341,22 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): material_id="pec", ) b16 = gprMax.Box( - p1=(x + 0.032 + cavitysize[0], y + 0.022, z + skidthickness + cavitysize[2] - casethickness), + p1=( + x + 0.032 + cavitysize[0], + y + 0.022, + z + skidthickness + cavitysize[2] - casethickness, + ), p2=(x + 0.108, y + 0.022 + 0.006, z + skidthickness + cavitysize[2]), material_id="pec", ) # PCB - replace bits chopped by TX & Rx cavities b17 = gprMax.Box( - p1=(x + 0.032 + cavitythickness, y + 0.022 + cavitythickness, z + skidthickness), + p1=( + x + 0.032 + cavitythickness, + y + 0.022 + cavitythickness, + z + skidthickness, + ), p2=( x + 0.032 + cavitysize[0] - cavitythickness, y + 0.022 + cavitysize[1] - cavitythickness, @@ -280,7 +365,11 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): material_id="pcb", ) b18 = gprMax.Box( - p1=(x + 0.108 + cavitythickness, y + 0.022 + cavitythickness, z + skidthickness), + p1=( + x + 0.108 + cavitythickness, + y + 0.022 + cavitythickness, + z + skidthickness, + ), p2=( x + 0.108 + cavitysize[0] - cavitythickness, y + 0.022 + cavitysize[1] - cavitythickness, @@ -300,7 +389,11 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): thickness=0, material_id="pec", ) - e1 = gprMax.Edge(p1=(tx[0], tx[1] - 0.001, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec") + e1 = gprMax.Edge( + p1=(tx[0], tx[1] - 0.001, tx[2]), + p2=(tx[0], tx[1], tx[2]), + material_id="pec", + ) t2 = gprMax.Triangle( p1=(tx[0], tx[1] + 0.002, tx[2]), p2=(tx[0] - 0.026, tx[1] + bowtieheight + 0.002, tx[2]), @@ -308,7 +401,11 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): thickness=0, material_id="pec", ) - e2 = gprMax.Edge(p1=(tx[0], tx[1] + 0.001, tx[2]), p2=(tx[0], tx[1] + 0.002, tx[2]), material_id="pec") + e2 = gprMax.Edge( + p1=(tx[0], tx[1] + 0.001, tx[2]), + p2=(tx[0], tx[1] + 0.002, tx[2]), + material_id="pec", + ) scene_objects.extend((t1, t2, e1, e2)) elif resolution == 0.002: t1 = gprMax.Triangle( @@ -336,7 +433,11 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): thickness=0, material_id="pec", ) - e3 = gprMax.Edge(p1=(tx[0] + 0.076, tx[1] - 0.001, tx[2]), p2=(tx[0] + 0.076, tx[1], tx[2]), material_id="pec") + e3 = gprMax.Edge( + p1=(tx[0] + 0.076, tx[1] - 0.001, tx[2]), + p2=(tx[0] + 0.076, tx[1], tx[2]), + material_id="pec", + ) t4 = gprMax.Triangle( p1=(tx[0] + 0.076, tx[1] + 0.002, tx[2]), p2=(tx[0] + 0.076 - 0.026, tx[1] + bowtieheight + 0.002, tx[2]), @@ -345,7 +446,9 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): material_id="pec", ) e4 = gprMax.Edge( - p1=(tx[0] + 0.076, tx[1] + 0.001, tx[2]), p2=(tx[0] + 0.076, tx[1] + 0.002, tx[2]), material_id="pec" + p1=(tx[0] + 0.076, tx[1] + 0.001, tx[2]), + p2=(tx[0] + 0.076, tx[1] + 0.002, tx[2]), + material_id="pec", ) scene_objects.extend((t3, e3, t4, e4)) elif resolution == 0.002: @@ -631,11 +734,17 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): # Skid b19 = gprMax.Box( - p1=(x, y, z), p2=(x + casesize[0], y + casesize[1], z + polypropylenethickness), material_id="polypropylene" + p1=(x, y, z), + p2=(x + casesize[0], y + casesize[1], z + polypropylenethickness), + material_id="polypropylene", ) b20 = gprMax.Box( p1=(x, y, z + polypropylenethickness), - p2=(x + casesize[0], y + casesize[1], z + polypropylenethickness + hdpethickness), + p2=( + x + casesize[0], + y + casesize[1], + z + polypropylenethickness + hdpethickness, + ), material_id="hdpe", ) scene_objects.extend((b19, b20)) @@ -644,7 +753,10 @@ def antenna_like_MALA_1200(x, y, z, resolution=0.001, **kwargs): w2 = gprMax.Waveform(wave_type="gaussian", amp=1, freq=excitationfreq, id="my_gaussian") scene_objects.append(w2) vs1 = gprMax.VoltageSource( - polarisation="y", p1=(tx[0], tx[1], tx[2]), resistance=sourceresistance, waveform_id="my_gaussian" + polarisation="y", + p1=(tx[0], tx[1], tx[2]), + resistance=sourceresistance, + waveform_id="my_gaussian", ) scene_objects.append(vs1) @@ -655,7 +767,11 @@ def antenna_like_MALA_1200(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_MALA_1200", output_type="n", diff --git a/toolboxes/Plotting/plot_Ascan.py b/toolboxes/Plotting/plot_Ascan.py index 25ad1c7c..f6db4bfe 100644 --- a/toolboxes/Plotting/plot_Ascan.py +++ b/toolboxes/Plotting/plot_Ascan.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. @@ -173,7 +173,10 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False, save=False): # Plotting if no FFT required else: fig, ax = plt.subplots( - subplot_kw=dict(xlabel="Time [s]", ylabel=outputtext + " field strength [V/m]"), + subplot_kw=dict( + xlabel="Time [s]", + ylabel=outputtext + " field strength [V/m]", + ), num=rxpath + " - " + f[rxpath].attrs["Name"], figsize=(20, 10), facecolor="w", @@ -279,7 +282,13 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False, save=False): if save: # Save a PDF of the figure - fig.savefig(filename[:-3] + ".pdf", dpi=None, format="pdf", bbox_inches="tight", pad_inches=0.1) + fig.savefig( + filename[:-3] + ".pdf", + dpi=None, + format="pdf", + bbox_inches="tight", + pad_inches=0.1, + ) # Save a PNG of the figure # fig.savefig(filename[:-3] + '.png', dpi=150, format='png', # bbox_inches='tight', pad_inches=0.1) @@ -322,9 +331,17 @@ if __name__ == "__main__": ], nargs="+", ) - parser.add_argument("-fft", action="store_true", default=False, help="plot FFT (single output must be specified)") parser.add_argument( - "-save", action="store_true", default=False, help="save plot directly to file, i.e. do not display" + "-fft", + action="store_true", + default=False, + help="plot FFT (single output must be specified)", + ) + parser.add_argument( + "-save", + action="store_true", + default=False, + help="save plot directly to file, i.e. do not display", ) args = parser.parse_args() diff --git a/toolboxes/Plotting/plot_Bscan.py b/toolboxes/Plotting/plot_Bscan.py index a2abb906..a42c482c 100644 --- a/toolboxes/Plotting/plot_Bscan.py +++ b/toolboxes/Plotting/plot_Bscan.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. @@ -46,7 +46,12 @@ def mpl_plot(filename, outputdata, dt, rxnumber, rxcomponent, save=False): file = Path(filename) - fig = plt.figure(num=file.stem + " - rx" + str(rxnumber), figsize=(20, 10), facecolor="w", edgecolor="w") + fig = plt.figure( + num=file.stem + " - rx" + str(rxnumber), + figsize=(20, 10), + facecolor="w", + edgecolor="w", + ) plt.imshow( outputdata, extent=[0, outputdata.shape[1], outputdata.shape[0] * dt, 0], @@ -73,7 +78,13 @@ def mpl_plot(filename, outputdata, dt, rxnumber, rxcomponent, save=False): if save: # Save a PDF of the figure - fig.savefig(filename[:-3] + ".pdf", dpi=None, format="pdf", bbox_inches="tight", pad_inches=0.1) + fig.savefig( + filename[:-3] + ".pdf", + dpi=None, + format="pdf", + bbox_inches="tight", + pad_inches=0.1, + ) # Save a PNG of the figure # fig.savefig(filename[:-3] + '.png', dpi=150, format='png', # bbox_inches='tight', pad_inches=0.1) @@ -94,10 +105,16 @@ if __name__ == "__main__": choices=["Ex", "Ey", "Ez", "Hx", "Hy", "Hz", "Ix", "Iy", "Iz"], ) parser.add_argument( - "-gather", action="store_true", default=False, help="gather together all receiver outputs in file" + "-gather", + action="store_true", + default=False, + help="gather together all receiver outputs in file", ) parser.add_argument( - "-save", action="store_true", default=False, help="save plot directly to file, i.e. do not display" + "-save", + action="store_true", + default=False, + help="save plot directly to file, i.e. do not display", ) args = parser.parse_args() @@ -118,7 +135,9 @@ if __name__ == "__main__": rxsgather = outputdata rxsgather = np.column_stack((rxsgather, outputdata)) else: - plthandle = mpl_plot(args.outputfile, outputdata, dt, rx, args.rx_component, save=args.save) + plthandle = mpl_plot( + args.outputfile, outputdata, dt, rx, args.rx_component, save=args.save + ) # Plot all receivers from single output file together if required if args.gather: diff --git a/toolboxes/Plotting/plot_antenna_params.py b/toolboxes/Plotting/plot_antenna_params.py index 08553cf2..6d7b9745 100644 --- a/toolboxes/Plotting/plot_antenna_params.py +++ b/toolboxes/Plotting/plot_antenna_params.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. diff --git a/toolboxes/Plotting/plot_source_wave.py b/toolboxes/Plotting/plot_source_wave.py index af89d153..c244a293 100644 --- a/toolboxes/Plotting/plot_source_wave.py +++ b/toolboxes/Plotting/plot_source_wave.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. @@ -90,7 +90,13 @@ def mpl_plot(w, timewindow, dt, iterations, fft=False, save=False): if w.freq and w.type != "gaussian" and w.type != "impulse": logging.info(f"Centre frequency: {w.freq:g} Hz") - if w.type in ["gaussian", "gaussiandot", "gaussiandotnorm", "gaussianprime", "gaussiandoubleprime"]: + if w.type in [ + "gaussian", + "gaussiandot", + "gaussiandotnorm", + "gaussianprime", + "gaussiandoubleprime", + ]: delay = 1 / w.freq logging.info(f"Time to centre of pulse: {delay:g} s") elif w.type in ["gaussiandotdot", "gaussiandotdotnorm", "ricker"]: @@ -113,7 +119,9 @@ def mpl_plot(w, timewindow, dt, iterations, fft=False, save=False): pltrange = np.where(freqs > 4 * w.freq)[0][0] pltrange = np.s_[0:pltrange] - fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, num=w.type, figsize=(20, 10), facecolor="w", edgecolor="w") + fig, (ax1, ax2) = plt.subplots( + nrows=1, ncols=2, num=w.type, figsize=(20, 10), facecolor="w", edgecolor="w" + ) # Plot waveform ax1.plot(time, waveform, "r", lw=2) @@ -143,9 +151,21 @@ def mpl_plot(w, timewindow, dt, iterations, fft=False, save=False): if save: savefile = Path(__file__).parent / w.type # Save a PDF of the figure - fig.savefig(savefile.with_suffix(".pdf"), dpi=None, format="pdf", bbox_inches="tight", pad_inches=0.1) + fig.savefig( + savefile.with_suffix(".pdf"), + dpi=None, + format="pdf", + bbox_inches="tight", + pad_inches=0.1, + ) # Save a PNG of the figure - fig.savefig(savefile.with_suffix(".png"), dpi=150, format="png", bbox_inches="tight", pad_inches=0.1) + fig.savefig( + savefile.with_suffix(".png"), + dpi=150, + format="png", + bbox_inches="tight", + pad_inches=0.1, + ) return plt @@ -163,16 +183,23 @@ if __name__ == "__main__": parser.add_argument("dt", type=float, help="time step to view waveform") parser.add_argument("-fft", action="store_true", default=False, help="plot FFT of waveform") parser.add_argument( - "-save", action="store_true", default=False, help="save plot directly to file, i.e. do not display" + "-save", + action="store_true", + default=False, + help="save plot directly to file, i.e. do not display", ) args = parser.parse_args() # Check waveform parameters if args.type.lower() not in Waveform.types: - logging.exception(f"The waveform must have one of the following types {', '.join(Waveform.types)}") + logging.exception( + f"The waveform must have one of the following types {', '.join(Waveform.types)}" + ) raise ValueError if args.freq <= 0: - logging.exception("The waveform requires an excitation frequency value of greater than zero") + logging.exception( + "The waveform requires an excitation frequency value of greater than zero" + ) raise ValueError # Create waveform instance diff --git a/toolboxes/STLtoVoxel/convert.py b/toolboxes/STLtoVoxel/convert.py index dfe02736..f3a383d3 100644 --- a/toolboxes/STLtoVoxel/convert.py +++ b/toolboxes/STLtoVoxel/convert.py @@ -26,7 +26,13 @@ def convert_files(input_file_paths, discretization, colors=[(0, 0, 0)], pad=1, p for input_file_path in input_file_paths: mesh_obj = mesh.Mesh.from_file(input_file_path) - org_mesh = np.hstack((mesh_obj.v0[:, np.newaxis], mesh_obj.v1[:, np.newaxis], mesh_obj.v2[:, np.newaxis])) + org_mesh = np.hstack( + ( + mesh_obj.v0[:, np.newaxis], + mesh_obj.v1[:, np.newaxis], + mesh_obj.v2[:, np.newaxis], + ) + ) meshes.append(org_mesh) vol, scale, shift = convert_meshes(meshes, discretization, parallel) vol = np.transpose(vol) diff --git a/toolboxes/STLtoVoxel/slice.py b/toolboxes/STLtoVoxel/slice.py index 96460d02..dc182527 100644 --- a/toolboxes/STLtoVoxel/slice.py +++ b/toolboxes/STLtoVoxel/slice.py @@ -19,7 +19,12 @@ def mesh_to_plane(mesh, bounding_box, parallel): current_mesh_indices = set() z = 0 - with tqdm(total=bounding_box[2], desc="Processing Layers", ncols=get_terminal_width() - 1, file=sys.stdout) as pbar: + with tqdm( + total=bounding_box[2], + desc="Processing Layers", + ncols=get_terminal_width() - 1, + file=sys.stdout, + ) as pbar: for event_z, status, tri_ind in generate_tri_events(mesh): while event_z - z >= 0: mesh_subset = [mesh[ind] for ind in current_mesh_indices] diff --git a/toolboxes/STLtoVoxel/stltovoxel.py b/toolboxes/STLtoVoxel/stltovoxel.py index 258eb970..bc323b98 100644 --- a/toolboxes/STLtoVoxel/stltovoxel.py +++ b/toolboxes/STLtoVoxel/stltovoxel.py @@ -17,9 +17,15 @@ if __name__ == "__main__": usage="cd gprMax; python -m toolboxes.STLtoVoxel.stltovoxel stlfilename -matindex -dxdydz", ) parser.add_argument( - "stlfiles", help="can be the filename of a single STL file, or the path to folder containing multiple STL files" + "stlfiles", + help="can be the filename of a single STL file, or the path to folder containing multiple STL files", + ) + parser.add_argument( + "-dxdydz", + type=float, + required=True, + help="discretisation to use in voxelisation process", ) - parser.add_argument("-dxdydz", type=float, required=True, help="discretisation to use in voxelisation process") args = parser.parse_args() if os.path.isdir(args.stlfiles): @@ -29,16 +35,18 @@ if __name__ == "__main__": filename_mats = os.path.join(path, os.path.basename(path) + "_mats.txt") elif os.path.isfile(args.stlfiles): path = os.path.dirname(args.stlfiles) - files = args.stlfiles - filename_hdf5 = os.path.join(path, os.path.split(os.path.basename(path))[0] + "_geo.h5") - filename_mats = os.path.join(path, os.path.split(os.path.basename(path))[0] + "_mats.txt") + files = [args.stlfiles] + filename_hdf5 = os.path.splitext(args.stlfiles)[0] + "_geo.h5" + filename_mats = os.path.splitext(args.stlfiles)[0] + "_mats.txt" dxdydz = (args.dxdydz, args.dxdydz, args.dxdydz) newline = "\n\t" logger.info(f"\nConverting STL file(s): {newline.join(files)}") model_array = convert_files(files, dxdydz) - logger.info(f"Number of voxels: {model_array.shape[0]} x {model_array.shape[1]} x {model_array.shape[2]}") + logger.info( + f"Number of voxels: {model_array.shape[0]} x {model_array.shape[1]} x {model_array.shape[2]}" + ) logger.info(f"Spatial discretisation: {dxdydz[0]} x {dxdydz[1]} x {dxdydz[2]}m") # Write HDF5 file for gprMax using voxels diff --git a/toolboxes/Utilities/Paraview/gprMax.py b/toolboxes/Utilities/Paraview/gprMax.py index 3adaf67a..4eafc7a3 100644 --- a/toolboxes/Utilities/Paraview/gprMax.py +++ b/toolboxes/Utilities/Paraview/gprMax.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. diff --git a/toolboxes/Utilities/convert_png2h5.py b/toolboxes/Utilities/convert_png2h5.py index 333220f9..32c84608 100644 --- a/toolboxes/Utilities/convert_png2h5.py +++ b/toolboxes/Utilities/convert_png2h5.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. @@ -55,7 +55,9 @@ class Cursor(object): ) # Convert pixel values from float (0-1) to integer (0-255) match = pixel_match(materials, pixel) if match is False: - logger.info(f"x, y: {int(x)} {int(y)} px; RGB: {pixel[:-1]}; material ID: {len(self.materials)}") + logger.info( + f"x, y: {int(x)} {int(y)} px; RGB: {pixel[:-1]}; material ID: {len(self.materials)}" + ) materials.append(pixel) @@ -85,10 +87,17 @@ if __name__ == "__main__": ) parser.add_argument("imagefile", help="name of image file including path") parser.add_argument( - "dxdydz", type=float, action="append", nargs=3, help="spatial resolution of model, e.g. dx dy dz" + "dxdydz", + type=float, + action="append", + nargs=3, + help="spatial resolution of model, e.g. dx dy dz", ) parser.add_argument( - "-zcells", default=1, type=int, help="number of cells for domain in z-direction (infinite direction)" + "-zcells", + default=1, + type=int, + help="number of cells for domain in z-direction (infinite direction)", ) args = parser.parse_args() @@ -97,7 +106,9 @@ if __name__ == "__main__": # Store image data to use for creating geometry imdata = np.rot90(im, k=3) # Rotate 90CW - imdata = np.floor(imdata * 255).astype(np.int16) # Convert pixel values from float (0-1) to integer (0-255) + imdata = np.floor(imdata * 255).astype( + np.int16 + ) # Convert pixel values from float (0-1) to integer (0-255) logger.info(f"Reading PNG image file: {os.path.split(args.imagefile)[1]}") logger.info( diff --git a/toolboxes/Utilities/get_host_spec.py b/toolboxes/Utilities/get_host_spec.py index 184dfa3c..ca392ce2 100644 --- a/toolboxes/Utilities/get_host_spec.py +++ b/toolboxes/Utilities/get_host_spec.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 @@ logging.basicConfig(format="%(message)s", level=logging.INFO) # Host machine info. hostinfo = get_host_info() -hyperthreadingstr = f", {hostinfo['logicalcores']} cores with Hyper-Threading" if hostinfo["hyperthreading"] else "" +hyperthreadingstr = ( + f", {hostinfo['logicalcores']} cores with Hyper-Threading" if hostinfo["hyperthreading"] else "" +) hostname = f"\n=== {hostinfo['hostname']}" logging.info(f"{hostname} {'=' * (get_terminal_width() - len(hostname) - 1)}") logging.info(f"\n{'Mfr/model:':<12} {hostinfo['machineID']}") @@ -46,7 +48,8 @@ logging.info(f"{'OS/Version:':<12} {hostinfo['osversion']}") # OpenMP logging.info( - "\n\n=== OpenMP capabilities (gprMax will not use Hyper-Threading " + "as there is no performance advantage)\n" + "\n\n=== OpenMP capabilities (gprMax will not use Hyper-Threading " + + "as there is no performance advantage)\n" ) logging.info(f"{'OpenMP threads: '} {hostinfo['physicalcores']}") diff --git a/toolboxes/Utilities/outputfiles_merge.py b/toolboxes/Utilities/outputfiles_merge.py index 2887f0a6..0b30539d 100644 --- a/toolboxes/Utilities/outputfiles_merge.py +++ b/toolboxes/Utilities/outputfiles_merge.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. @@ -110,7 +110,9 @@ def merge_files(outputfiles, merged_outputfile=None, removefiles=False): availableoutputs = list(fin[path].keys()) for output in availableoutputs: grp.create_dataset( - output, (fout.attrs["Iterations"], len(outputfiles)), dtype=fin[path + "/" + output].dtype + output, + (fout.attrs["Iterations"], len(outputfiles)), + dtype=fin[path + "/" + output].dtype, ) # For all receivers @@ -139,14 +141,24 @@ if __name__ == "__main__": ) parser.add_argument("basefilename", help="base name of output file series including path") parser.add_argument( - "-o", "--output-file", default=None, type=str, required=False, help="location to save merged file" + "-o", + "--output-file", + default=None, + type=str, + required=False, + help="location to save merged file", ) parser.add_argument( - "--remove-files", action="store_true", default=False, help="flag to remove individual output files after merge" + "--remove-files", + action="store_true", + default=False, + help="flag to remove individual output files after merge", ) args = parser.parse_args() files = glob.glob(args.basefilename + "*.h5") - outputfiles = [filename for filename in files if "_merged" not in filename and args.output_file != filename] + outputfiles = [ + filename for filename in files if "_merged" not in filename and args.output_file != filename + ] outputfiles.sort(key=natural_keys) merge_files(outputfiles, merged_outputfile=args.output_file, removefiles=args.remove_files)