diff --git a/README.rst b/README.rst index 07b3c31b..8b39624d 100644 --- a/README.rst +++ b/README.rst @@ -117,7 +117,7 @@ Microsoft Windows ^^^^^^^^^^^^^^^^^ * Download and install Microsoft `Build Tools for Visual Studio 2022 `_ (direct link). You can also find it on the `Microsoft Visual Studio downloads page `_ by scrolling down to the 'All Downloads' section, clicking the disclosure triangle by 'Tools for Visual Studio 2022', then clicking the download button next to 'Build Tools for Visual Studio 2022'. When installing, choose the 'Desktop development with C++' Workload and select only 'MSVC v143' and 'Windows 10 SDK' or 'Windows 11 SDK options. -* Set the Path and Environment Variables - this can be done by following the `instructions from Microsoft `_, or manually by adding a form of :code:`C:\Program Files (x86)\Microsoft Visual Studio\2019\BuildTools\VC\Tools\MSVC\14.23.28105\bin\Hostx64\x64` (this may vary according to your exact machine and installation) to your system Path environment variable. +* Set the Path and Environment Variables - this can be done by following the `instructions from Microsoft `_, or manually by adding a form of :code:``C:\Program Files (x86)\Microsoft Visual Studio\2019\BuildTools\VC\Tools\MSVC\14.23.28105\bin\Hostx64\x64`` (this may vary according to your exact machine and installation) to your system Path environment variable. Alternatively, if you are using Windows 10/11 you can install the `Windows Subsystem for Linux `_ and then follow the Linux install instructions for gprMax. Note however that currently, WSL does not aim to support GUI desktops or applications, e.g. Gnome, KDE, etc.... diff --git a/docs/source/examples_simple_2D.rst b/docs/source/examples_simple_2D.rst index 0ed67e0d..c42ed962 100644 --- a/docs/source/examples_simple_2D.rst +++ b/docs/source/examples_simple_2D.rst @@ -112,7 +112,7 @@ You should have produced an output file ``cylinder_Ascan_2D.h5``. You can view t .. code-block:: none - python -m tools.plot_Ascan examples/cylinder_Ascan_2D.h5 + python -m toolboxes.Plotting.plot_Ascan examples/cylinder_Ascan_2D.h5 :numref:`cylinder_Ascan_results` shows the time history of the electric and magnetic field components and currents at the receiver location. The :math:`E_z` field component can be converted to a voltage that represents the A-scan (trace). The initial part of the signal (~0.5-1.5 ns) represents the direct wave from transmitter to receiver. Then comes the reflected wavelet (~1.8-2.6 ns), which has opposite polarity, from the metal cylinder. @@ -153,7 +153,7 @@ You should have produced 60 output files, one for each A-scan, with names ``cyli .. code-block:: none - python -m tools.outputfiles_merge examples/cylinder_Bscan_2D + python -m toolboxes.Utilities.outputfiles_merge examples/cylinder_Bscan_2D You should see a combined output file ``cylinder_Bscan_2D_merged.h5``. You can add the optional argument ``--remove-files`` if you want to automatically delete the original single A-scan output files. @@ -161,7 +161,7 @@ You can now view an image of the B-scan using the command: .. code-block:: none - python -m tools.plot_Bscan examples/cylinder_Bscan_2D_merged.h5 Ez + python -m toolboxes.Plotting.plot_Bscan examples/cylinder_Bscan_2D_merged.h5 Ez :numref:`cylinder_Bscan_results` shows the B-scan (of the :math:`E_z` field component). Again, the initial part of the signal (~0.5-1.5 ns) represents the direct wave from transmitter to receiver. Then comes the reflected wave (~2-3 ns) from the metal cylinder which creates the hyperbolic shape. diff --git a/docs/source/input_hash_cmds.rst b/docs/source/input_hash_cmds.rst index c4710f32..e6ebf26a 100644 --- a/docs/source/input_hash_cmds.rst +++ b/docs/source/input_hash_cmds.rst @@ -965,7 +965,7 @@ For example to save a snapshot of the electromagnetic fields in the model at a s PML commands ============ -The default behaviour for the absorbing boundary conditions (ABC) is first order Complex Frequency Shifted (CFS) Perfectly Matched Layers (PML), with thicknesses of 10 cells on each of the six sides of the model domain. The thickness of the PML can be altered by using the following command (further customisation of the PML is possible using our Python API): +The default behaviour for the absorbing boundary conditions (ABC) is first order Complex Frequency Shifted (CFS) Perfectly Matched Layers (PML), with thicknesses of 10 cells on each of the six sides of the model domain. #pml_cells: ------------ @@ -989,3 +989,47 @@ For example to use a PML with 20 cells (thicker than the default 10 cells) on on .. code-block:: none #pml_cells: 10 10 20 10 10 20 + +#pml_formulation: +----------------- + +Allows you to alter the formulation used for the PML. The current options are to use the Higher Order RIPML (HORIPML) - https://doi.org/10.1109/TAP.2011.2180344, or Multipole RIPML (MRIPML) - https://doi.org/10.1109/TAP.2018.2823864. The syntax of the command is: + +.. code-block:: none + + #pml_formulation: str + +* ``str`` can be either 'HORIPML' or 'MRIPML' + +For example to use the Multipole RIPML: + +.. code-block:: none + + #pml_formulation: MRIPML + +#pml_cfs: +--------- + +Allows you (advanced) control of the parameters that are used to build each order of the PML. Up to a second order PML can currently be specified, i.e. by using two ``#pml_cfs`` commands. The syntax of the command is: + +.. code-block:: none + + #pml_cfs: str1 str2 f1 f2 str3 str4 f3 f4 str5 str6 f5 f6 + +* ``str1`` is the type of scaling to use for the CFS :math:`\alpha` parameter. It can be ``constant``, ``linear``, ``quadratic``, ``cubic``, ``quartic``, ``quintic`` and ``sextic``. +* ``str2`` is the direction of the scaling to use for the CFS :math:`\alpha` parameter. It can be ``forward`` or ``reverse``. +* ``f1 f2`` are the minimum and maximum values for the CFS :math:`\alpha` parameter. +* ``str3`` is the type of scaling to use for the CFS :math:`\kappa` parameter. It can be ``constant``, ``linear``, ``quadratic``, ``cubic``, ``quartic``, ``quintic`` and ``sextic``. +* ``str4`` is the direction of the scaling to use for the CFS :math:`\kappa` parameter. It can be ``forward`` or ``reverse``. +* ``f3 f4`` are the minimum and maximum values for the CFS :math:`\kappa` parameter. The minimum value for the CFS :math:`\kappa` parameter is one. +* ``str5`` is the type of scaling to use for the CFS :math:`\sigma` parameter. It can be ``constant``, ``linear``, ``quadratic``, ``cubic``, ``quartic``, ``quintic`` and ``sextic``. +* ``str6`` is the direction of the scaling to use for the CFS :math:`\sigma` parameter. It can be ``forward`` or ``reverse``. +* ``f5 f6`` are the minimum and maximum values for the CFS :math:`\sigma` parameter. + +The CFS values (which are internally specified) used for the default standard first order PML are: ``#pml_cfs: constant forward 0 0 constant forward 1 1 quartic forward 0 None``. Specifying 'None' for the maximum value of :math:`\sigma` forces gprMax to calculate it internally based on the relative permittivity and permeability of the underlying materials in the model. + +The parameters will be applied to all slabs of the PML that are switched on. + +.. tip:: + + ``forward`` direction implies minimum parameter value at the inner boundary of the PML and maximum parameter value at the edge of computational domain, ``reverse`` is the opposite. diff --git a/examples/antenna_like_GSSI_400_fs.py b/examples/antenna_like_GSSI_400_fs.py index 1477c12f..56c52b6b 100644 --- a/examples/antenna_like_GSSI_400_fs.py +++ b/examples/antenna_like_GSSI_400_fs.py @@ -30,7 +30,7 @@ scene = gprMax.Scene() title = gprMax.Title(name=fn.with_suffix("").name) domain = gprMax.Domain(p1=(x, y, z)) dxdydz = gprMax.Discretisation(p1=(dl, dl, dl)) -time_window = gprMax.TimeWindow(time=6e-9) +time_window = gprMax.TimeWindow(time=15e-9) scene.add(title) scene.add(domain) @@ -44,17 +44,20 @@ for obj in gssi_objects: scene.add(obj) gv1 = gprMax.GeometryView( - p1=(0, 0, 0), p2=(x, y, z), dl=(dl, dl, dl), filename="antenna_like_GSSI_400", output_type="n" + p1=(0, 0, 0), + p2=(x, y, z), + dl=(dl, dl, dl), + filename="antenna_like_GSSI_400", output_type="n" ) gv2 = gprMax.GeometryView( - p1=(ant_pos[0] - 0.150 / 2, ant_pos[1] - 0.150 / 2, ant_pos[2] - 0.050), - p2=(ant_pos[0] + 0.150 / 2, ant_pos[1] + 0.150 / 2, ant_pos[2] + 0.010), + p1=(ant_pos[0] - 0.150, ant_pos[1] - 0.150, ant_pos[2]), + p2=(ant_pos[0] + 0.150, ant_pos[1] + 0.150, ant_pos[2] + 0.050), dl=(dl, dl, dl), filename="antenna_like_GSSI_400_pcb", output_type="f", ) -scene.add(gv1) -scene.add(gv2) +#scene.add(gv1) +#scene.add(gv2) # Run model gprMax.run(scenes=[scene], geometry_only=False, outputfile=fn, gpu=None) diff --git a/gprMax/__init__.py b/gprMax/__init__.py index 7c48639f..4e572e47 100644 --- a/gprMax/__init__.py +++ b/gprMax/__init__.py @@ -14,9 +14,11 @@ from .cmds_geometry.add_grass import AddGrass from .cmds_geometry.add_surface_roughness import AddSurfaceRoughness from .cmds_geometry.add_surface_water import AddSurfaceWater from .cmds_geometry.box import Box +from .cmds_geometry.cone import Cone from .cmds_geometry.cylinder import Cylinder from .cmds_geometry.cylindrical_sector import CylindricalSector from .cmds_geometry.edge import Edge +from .cmds_geometry.ellipsoid import Ellipsoid from .cmds_geometry.fractal_box import FractalBox from .cmds_geometry.geometry_objects_read import GeometryObjectsRead from .cmds_geometry.plate import Plate @@ -38,6 +40,7 @@ from .cmds_multiuse import ( RxArray, Snapshot, SoilPeplinski, + Subgrid, TransmissionLine, VoltageSource, Waveform, @@ -47,6 +50,7 @@ from .cmds_singleuse import ( Domain, ExcitationFile, OMPThreads, + OutputDir, PMLProps, RxSteps, SrcSteps, diff --git a/gprMax/cmds_multiuse.py b/gprMax/cmds_multiuse.py index 378d3261..0388c74f 100644 --- a/gprMax/cmds_multiuse.py +++ b/gprMax/cmds_multiuse.py @@ -1737,11 +1737,6 @@ class PMLCFS(UserObjectMulti): ): logger.exception(f"{self.params_str()} minimum and maximum scaling values must be greater than zero.") raise ValueError - # TODO: Fix handling of kappa for 2nd order PMLs - # if float(kappamin) < 1: - # logger.exception(f'{self.params_str()} minimum scaling value for ' - # 'kappa must be greater than or equal to one.') - # raise ValueError cfsalpha = CFSParameter() cfsalpha.ID = "alpha" @@ -1760,8 +1755,11 @@ class PMLCFS(UserObjectMulti): cfssigma.scalingprofile = sigmascalingprofile cfssigma.scalingdirection = sigmascalingdirection cfssigma.min = float(sigmamin) + if sigmamax == "None": + sigmamax = None if sigmamax is not None: - cfssigma.max = float(sigmamax) + sigmamax = float(sigmamax) + cfssigma.max = sigmamax cfs = CFS() cfs.alpha = cfsalpha cfs.kappa = cfskappa @@ -1775,7 +1773,7 @@ class PMLCFS(UserObjectMulti): f"{cfskappa.scalingdirection}, min: {cfskappa.min:g}, max: " f"{cfskappa.max:g}), sigma (scaling: {cfssigma.scalingprofile}, " f"scaling direction: {cfssigma.scalingdirection}, min: " - f"{cfssigma.min:g}, max: {cfssigma.max:g}) created." + f"{cfssigma.min:g}, max: {cfssigma.max}) created." ) grid.pmls["cfs"].append(cfs) @@ -1800,4 +1798,4 @@ class Subgrid(UserObjectMulti): self.children_geometry.append(node) else: logger.exception("This object is unknown to gprMax.") - raise ValueError + raise ValueError diff --git a/gprMax/cmds_singleuse.py b/gprMax/cmds_singleuse.py index 2570bb86..b2560548 100644 --- a/gprMax/cmds_singleuse.py +++ b/gprMax/cmds_singleuse.py @@ -433,7 +433,7 @@ class ExcitationFile(UserObjectSingle): waveformtime = np.arange(0, G.timewindow + G.dt, G.dt) timestr = "simulation time array" - for waveformID in waveformIDs: + for i, waveformID in enumerate(waveformIDs): if any(x.ID == waveformID for x in G.waveforms): logger.exception(f"Waveform with ID {waveformID} already exists") raise ValueError @@ -442,7 +442,7 @@ class ExcitationFile(UserObjectSingle): w.type = "user" # Select correct column of waveform values depending on array shape - singlewaveformvalues = waveformvalues[:] if len(waveformvalues.shape) == 1 else waveformvalues[:, waveform] + singlewaveformvalues = waveformvalues[:] if len(waveformvalues.shape) == 1 else waveformvalues[:, i] # Truncate waveform array if it is longer than time array if len(singlewaveformvalues) > len(waveformtime): diff --git a/gprMax/config.py b/gprMax/config.py index fbbb7628..950072aa 100644 --- a/gprMax/config.py +++ b/gprMax/config.py @@ -209,33 +209,29 @@ class SimulationConfig: self.args = args - if args.mpi and args.geometry_fixed: + if self.args.mpi and self.args.geometry_fixed: logger.exception("The geometry fixed option cannot be used with MPI.") raise ValueError - if args.gpu and args.opencl: + if self.args.gpu and self.args.opencl: logger.exception("You cannot use both CUDA and OpenCl simultaneously.") raise ValueError # General settings for the simulation - # inputfilepath: path to inputfile location. - # outputfilepath: path to outputfile location. - # progressbars: whether to show progress bars on stdoout or not. # solver: cpu, cuda, opencl. - # subgrid: whether the simulation uses sub-grids. # precision: data type for electromagnetic field output (single/double). # progressbars: progress bars on stdoout or not - switch off # progressbars when logging level is greater than # info (20) - self.general = {"solver": "cpu", "subgrid": False, "precision": "single", "progressbars": args.log_level <= 20} + self.general = {"solver": "cpu", "precision": "single", "progressbars": args.log_level <= 20} self.em_consts = { "c": c, # Speed of light in free space (m/s) "e0": e0, # Permittivity of free space (F/m) "m0": m0, # Permeability of free space (H/m) - "z0": np.sqrt(m0 / e0), - } # Impedance of free space (Ohms) + "z0": np.sqrt(m0 / e0), # Impedance of free space (Ohms) + } # Store information about host machine self.hostinfo = get_host_info() @@ -246,7 +242,7 @@ class SimulationConfig: # Both single and double precision are possible on GPUs, but single # provides best performance. self.general["precision"] = "single" - self.devices = {"devs": [], "nvcc_opts": None} # pycuda device objects # nvcc compiler options + self.devices = {"devs": [], "nvcc_opts": None} # pycuda device objects; nvcc compiler options # Suppress nvcc warnings on Microsoft Windows if sys.platform == "win32": self.devices["nvcc_opts"] = ["-w"] @@ -258,7 +254,7 @@ class SimulationConfig: if self.args.opencl is not None: self.general["solver"] = "opencl" self.general["precision"] = "single" - self.devices = {"devs": [], "compiler_opts": None} # pyopencl available device(s) + self.devices = {"devs": [], "compiler_opts": None} # pyopencl device device(s); compiler options # Suppress unused variable warnings on gcc # if sys.platform != 'win32': self.devices['compiler_opts'] = ['-w'] @@ -266,8 +262,8 @@ class SimulationConfig: # Add pyopencl available device(s) self.devices["devs"] = detect_opencl() - # Subgrid parameter may not exist if user enters via CLI - try: + # Subgrids + if hasattr(self.args, "subgrid") and self.args.subgrid: self.general["subgrid"] = self.args.subgrid # Double precision should be used with subgrid for best accuracy self.general["precision"] = "double" @@ -278,7 +274,7 @@ class SimulationConfig: "You cannot currently use CUDA or OpenCL-based " "solvers with models that contain sub-grids." ) raise ValueError - except AttributeError: + else: self.general["subgrid"] = False # Scenes parameter may not exist if user enters via CLI diff --git a/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py b/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py index 8974a3a6..8d6a357e 100644 --- a/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py +++ b/gprMax/cuda_opencl/knl_pml_updates_electric_HORIPML.py @@ -35,11 +35,19 @@ x_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, $REAL *Ey, + const $REAL* __restrict__ Ex, + $REAL *Ey, $REAL *Ez, - const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, + const $REAL* __restrict__ Hx, + const $REAL* __restrict__ Hy, + const $REAL* __restrict__ Hz, + $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -95,9 +103,16 @@ y_args = { $REAL *Ex, const $REAL* __restrict__ Ey, $REAL *Ez, - const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, + const $REAL* __restrict__ Hx, + const $REAL* __restrict__ Hy, + const $REAL* __restrict__ Hz, + $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -153,9 +168,16 @@ z_args = { $REAL *Ex, $REAL *Ey, const $REAL* __restrict__ Ez, - const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, + const $REAL* __restrict__ Hx, + const $REAL* __restrict__ Hy, + const $REAL* __restrict__ Hz, + $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( diff --git a/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py b/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py index 1effe788..33adee2d 100644 --- a/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py +++ b/gprMax/cuda_opencl/knl_pml_updates_electric_MRIPML.py @@ -35,11 +35,19 @@ x_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, $REAL *Ey, + const $REAL* __restrict__ Ex, + $REAL *Ey, $REAL *Ez, - const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, + const $REAL* __restrict__ Hx, + const $REAL* __restrict__ Hy, + const $REAL* __restrict__ Hz, + $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -95,9 +103,16 @@ y_args = { $REAL *Ex, const $REAL* __restrict__ Ey, $REAL *Ez, - const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, + const $REAL* __restrict__ Hx, + const $REAL* __restrict__ Hy, + const $REAL* __restrict__ Hz, + $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -153,9 +168,16 @@ z_args = { $REAL *Ex, $REAL *Ey, const $REAL* __restrict__ Ez, - const $REAL* __restrict__ Hx, const $REAL* __restrict__ Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, + const $REAL* __restrict__ Hx, + const $REAL* __restrict__ Hy, + const $REAL* __restrict__ Hz, + $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( diff --git a/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py b/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py index 3f8bc823..6fccb4f6 100644 --- a/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py +++ b/gprMax/cuda_opencl/knl_pml_updates_magnetic_HORIPML.py @@ -35,11 +35,19 @@ x_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, const $REAL* __restrict__ Hx, $REAL *Hy, + const $REAL* __restrict__ Ex, + const $REAL* __restrict__ Ey, + const $REAL* __restrict__ Ez, + const $REAL* __restrict__ Hx, + $REAL *Hy, $REAL *Hz, $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -92,13 +100,19 @@ y_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, + const $REAL* __restrict__ Ex, + const $REAL* __restrict__ Ey, + const $REAL* __restrict__ Ez, $REAL *Hx, const $REAL* __restrict__ Hy, $REAL *Hz, $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -151,13 +165,19 @@ z_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, + const $REAL* __restrict__ Ex, + const $REAL* __restrict__ Ey, + const $REAL* __restrict__ Ez, $REAL *Hx, $REAL *Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( diff --git a/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py b/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py index 05ca8b9e..89df47ae 100644 --- a/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py +++ b/gprMax/cuda_opencl/knl_pml_updates_magnetic_MRIPML.py @@ -35,11 +35,19 @@ x_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, const $REAL* __restrict__ Hx, $REAL *Hy, + const $REAL* __restrict__ Ex, + const $REAL* __restrict__ Ey, + const $REAL* __restrict__ Ez, + const $REAL* __restrict__ Hx, + $REAL *Hy, $REAL *Hz, $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -92,13 +100,19 @@ y_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, + const $REAL* __restrict__ Ex, + const $REAL* __restrict__ Ey, + const $REAL* __restrict__ Ez, $REAL *Hx, const $REAL* __restrict__ Hy, $REAL *Hz, $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( @@ -151,13 +165,19 @@ z_args = { int NZ_PHI2, int NY_R, const unsigned int* __restrict__ ID, - const $REAL* __restrict__ Ex, const $REAL* __restrict__ Ey, const $REAL* __restrict__ Ez, + const $REAL* __restrict__ Ex, + const $REAL* __restrict__ Ey, + const $REAL* __restrict__ Ez, $REAL *Hx, $REAL *Hy, const $REAL* __restrict__ Hz, $REAL *PHI1, $REAL *PHI2, - const $REAL* __restrict__ RA, const $REAL* __restrict__ RB, const $REAL* __restrict__ RE, const $REAL* __restrict__ RF, $REAL d) + const $REAL* __restrict__ RA, + const $REAL* __restrict__ RB, + const $REAL* __restrict__ RE, + const $REAL* __restrict__ RF, + $REAL d) """ ), "opencl": Template( diff --git a/gprMax/cuda_opencl/knl_source_updates.py b/gprMax/cuda_opencl/knl_source_updates.py index 5e1e6372..85203c32 100644 --- a/gprMax/cuda_opencl/knl_source_updates.py +++ b/gprMax/cuda_opencl/knl_source_updates.py @@ -126,10 +126,10 @@ update_magnetic_dipole = { $REAL dx, $REAL dy, $REAL dz, - __global const int* restrict rcinfo1, - __global const $REAL* restrict rcinfo2, - __global const $REAL* restrict rcwaveforms, - __global const unsigned int* estrict ID, + __global const int* restrict srcinfo1, + __global const $REAL* restrict srcinfo2, + __global const $REAL* restrict srcwaveforms, + __global const unsigned int* restrict ID, __global $REAL *Hx, __global $REAL *Hy, __global $REAL *Hz @@ -208,10 +208,10 @@ update_voltage_source = { $REAL dx, $REAL dy, $REAL dz, - __global const int* restrict rcinfo1, - __global const $REAL* restrict rcinfo2, - __global const $REAL* restrict rcwaveforms, - __global const unsigned int* estrict ID, + __global const int* restrict srcinfo1, + __global const $REAL* restrict srcinfo2, + __global const $REAL* restrict srcwaveforms, + __global const unsigned int* restrict ID, __global $REAL *Ex, __global $REAL *Ey, __global $REAL *Ez diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index ba61add0..d3b31521 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.py @@ -25,7 +25,7 @@ from .utilities.logging import logging_config # Arguments (used for API) and their default values (used for API and CLI) args_defaults = { - "scenes": None, + "scenes": [], "inputfile": None, "outputfile": None, "n": 1, @@ -44,37 +44,39 @@ args_defaults = { # Argument help messages (used for CLI argparse) help_msg = { - "scenes": "(list, req): Scenes to run the model. " - "Multiple scene objects can given in order to run multiple " - "simulation runs. Each scene must contain the essential " - "simulation objects", - "inputfile": "(str, opt): Input file path. Can also run simulation " "by providing an input file.", + "scenes": "(list, req): Scenes to run the model. Multiple scene objects " + "can given in order to run multiple simulation runs. Each scene " + "must contain the essential simulation objects", + "inputfile": "(str, opt): Input file path. Can also run simulation by " + "providing an input file.", "outputfile": "(str, req): File path to the output data file.", "n": "(int, req): Number of required simulation runs.", - "i": "(int, opt): Model number to start/restart simulation " - "from. It would typically be used to restart a series of " - "models from a specific model number, with the n argument, " - "e.g. to restart from A-scan 45 when creating a B-scan " - "with 60 traces.", - "mpi": "(bool, opt): Flag to use Message Passing Interface (MPI) " - "task farm. This option is most usefully combined with n to " - "allow individual models to be farmed out using a MPI task " - "farm, e.g. to create a B-scan with 60 traces and use MPI to " - "farm out each trace. For further details see the parallel " - "performance section of the User Guide.", - "gpu": "(list/bool, opt): Flag to use NVIDIA GPU or list of NVIDIA " "GPU device ID(s) for specific GPU card(s).", - "opencl": "(list/bool, opt): Flag to use OpenCL or list of OpenCL " "device ID(s) for specific compute device(s).", + "i": "(int, opt): Model number to start/restart simulation from. It would " + "typically be used to restart a series of models from a specific " + "model number, with the n argument, e.g. to restart from A-scan 45 " + "when creating a B-scan with 60 traces.", + "mpi": "(bool, opt): Flag to use Message Passing Interface (MPI) task farm. " + "This option is most usefully combined with n to allow individual " + "models to be farmed out using a MPI task farm, e.g. to create a " + "B-scan with 60 traces and use MPI to farm out each trace. For " + "further details see the performance section of the User Guide.", + "gpu": "(list/bool, opt): Flag to use NVIDIA GPU or list of NVIDIA GPU " + "device ID(s) for specific GPU card(s).", + "opencl": "(list/bool, opt): Flag to use OpenCL or list of OpenCL device " + "ID(s) for specific compute device(s).", "subgrid": "(bool, opt): Flag to use sub-gridding.", - "autotranslate": "(bool, opt): For sub-gridding - auto translate " - "objects with main grid coordinates to their " - "equivalent local grid coordinate within the " - "subgrid. If this option is off users must specify " - "sub-grid object point within the global subgrid space.", - "geometry_only": "(bool, opt): Build a model and produce any " "geometry views but do not run the simulation.", - "geometry_fixed": "(bool, opt): Run a series of models where the " "geometry does not change between models.", - "write_processed": "(bool, opt): Writes another input file after " - "any Python code (#python blocks) and in the " - "original input file has been processed.", + "autotranslate": "(bool, opt): For sub-gridding - auto translate objects " + "with main grid coordinates to their equivalent local " + "grid coordinate within the subgrid. If this option is " + "off users must specify sub-grid object point within the " + "global subgrid space.", + "geometry_only": "(bool, opt): Build a model and produce any geometry " + "views but do not run the simulation.", + "geometry_fixed": "(bool, opt): Run a series of models where the geometry " + "does not change between models.", + "write_processed": "(bool, opt): Writes another input file after any " + "Python code (#python blocks) and in the original input " + "file has been processed.", "log_level": "(int, opt): Level of logging to use.", "log_file": "(bool, opt): Write logging information to file.", } @@ -209,9 +211,7 @@ def run_main(args): """ results = {} - logging_config(level=args.log_level, log_file=args.log_file) - config.sim_config = config.SimulationConfig(args) # MPI running with (OpenMP/CUDA/OpenCL) @@ -222,4 +222,5 @@ def run_main(args): context = Context() results = context.run() + return results diff --git a/gprMax/grid.py b/gprMax/grid.py index bb0ab22e..79141a45 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -56,10 +56,11 @@ class FDTDGrid: self.pmls["formulation"] = "HORIPML" self.pmls["cfs"] = [] self.pmls["slabs"] = [] - # Ordered dictionary required so that PMLs are always updated in the - # same order. The order itself does not matter, however, if must be the - # same from model to model otherwise the numerical precision from adding - # the PML corrections will be different. + # Ordered dictionary required so *updating* the PMLs always follows the + # same order (the order for *building* PMLs does not matter). The order + # itself does not matter, however, if must be the same from model to + # model otherwise the numerical precision from adding the PML + # corrections will be different. self.pmls["thickness"] = OrderedDict((key, 10) for key in PML.boundaryIDs) self.materials = [] @@ -353,6 +354,7 @@ class CUDAGrid(FDTDGrid): self.Hx_dev = gpuarray.to_gpu(self.Hx) self.Hy_dev = gpuarray.to_gpu(self.Hy) self.Hz_dev = gpuarray.to_gpu(self.Hz) + elif config.sim_config.general["solver"] == "opencl": import pyopencl.array as clarray diff --git a/gprMax/hash_cmds_file.py b/gprMax/hash_cmds_file.py index b5ba0bcb..90539a92 100644 --- a/gprMax/hash_cmds_file.py +++ b/gprMax/hash_cmds_file.py @@ -217,6 +217,7 @@ def check_cmd_names(processedlines, checkessential=True): "#title", "#omp_threads", "#time_step_stability_factor", + "#pml_formulation", "#pml_cells", "#excitation_file", "#src_steps", @@ -248,6 +249,7 @@ def check_cmd_names(processedlines, checkessential=True): "#rx", "#rx_array", "#snapshot", + "#pml_cfs", "#include_file", ] } diff --git a/gprMax/hash_cmds_multiuse.py b/gprMax/hash_cmds_multiuse.py index 18e71c06..f3e1f5ff 100644 --- a/gprMax/hash_cmds_multiuse.py +++ b/gprMax/hash_cmds_multiuse.py @@ -29,6 +29,7 @@ from .cmds_multiuse import ( Material, MaterialList, MaterialRange, + PMLCFS, Rx, RxArray, Snapshot, @@ -375,7 +376,7 @@ def process_multicmds(multicmds): tmp = cmdinstance.split() if len(tmp) < 2: - logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least 2 parameters") + logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires at least two parameters") raise ValueError tokens = len(tmp) @@ -386,4 +387,28 @@ def process_multicmds(multicmds): material_list = MaterialList(list_of_materials=lmats, id=tmp[tokens - 1]) scene_objects.append(material_list) + cmdname = '#pml_cfs' + if multicmds[cmdname] is not None: + for cmdinstance in multicmds[cmdname]: + tmp = cmdinstance.split() + + if len(tmp) != 12: + logger.exception("'" + cmdname + ": " + " ".join(tmp) + "'" + " requires exactly twelve parameters") + raise ValueError + + pml_cfs = PMLCFS(alphascalingprofile=tmp[0], + alphascalingdirection=tmp[1], + alphamin=tmp[2], + alphamax=tmp[3], + kappascalingprofile=tmp[4], + kappascalingdirection=tmp[5], + kappamin=tmp[6], + kappamax=tmp[7], + sigmascalingprofile=tmp[8], + sigmascalingdirection=tmp[9], + sigmamin=tmp[10], + sigmamax=tmp[11]) + + scene_objects.append(pml_cfs) + return scene_objects diff --git a/gprMax/hash_cmds_singleuse.py b/gprMax/hash_cmds_singleuse.py index 023e78a1..21fea96d 100644 --- a/gprMax/hash_cmds_singleuse.py +++ b/gprMax/hash_cmds_singleuse.py @@ -120,20 +120,40 @@ def process_singlecmds(singlecmds): scene_objects.append(tw) + cmd = "#pml_formulation" + if singlecmds[cmd] is not None: + tmp = singlecmds[cmd].split() + if len(tmp) != 1: + logger.exception(f"{cmd} requires one parameter") + raise ValueError + else: + pml_formulation = tmp[0] + cmd = "#pml_cells" if singlecmds[cmd] is not None: tmp = singlecmds[cmd].split() if len(tmp) not in [1, 6]: logger.exception(f"{cmd} requires either one or six parameter(s)") raise ValueError - if len(tmp) == 1: - pml_cells = PMLProps(thickness=int(tmp[0])) + + if 'pml_formulation' in locals(): + if len(tmp) == 1: + pml_props = PMLProps(formulation=pml_formulation, thickness=int(tmp[0])) + else: + pml_props = PMLProps( + formulation=pml_formulation, + x0=int(tmp[0]), y0=int(tmp[1]), z0=int(tmp[2]), + xmax=int(tmp[3]), ymax=int(tmp[4]), zmax=int(tmp[5]) + ) else: - pml_cells = PMLProps( - x0=int(tmp[0]), y0=int(tmp[1]), z0=int(tmp[2]), xmax=int(tmp[3]), ymax=int(tmp[4]), zmax=int(tmp[5]) - ) - - scene_objects.append(pml_cells) + if len(tmp) == 1: + pml_props = PMLProps(thickness=int(tmp[0])) + else: + pml_props = PMLProps( + x0=int(tmp[0]), y0=int(tmp[1]), z0=int(tmp[2]), xmax=int(tmp[3]), ymax=int(tmp[4]), zmax=int(tmp[5]) + ) + + scene_objects.append(pml_props) cmd = "#src_steps" if singlecmds[cmd] is not None: diff --git a/gprMax/pml.py b/gprMax/pml.py index fba68ddd..6ce24e5b 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.py @@ -433,19 +433,6 @@ class CUDAPML(PML): 1, ) - def get_update_funcs(self, kernelselectric, kernelsmagnetic): - """Gets update functions from PML kernels. - - Args: - kernelselectric: pycuda SourceModule containing PML kernels for - electric updates. - kernelsmagnetic: pycuda SourceModule containing PML kernels for - magnetic updates. - """ - - self.update_electric_dev = kernelselectric.get_function("order" + str(len(self.CFS)) + "_" + self.direction) - self.update_magnetic_dev = kernelsmagnetic.get_function("order" + str(len(self.CFS)) + "_" + self.direction) - def update_electric(self): """Updates electric field components with the PML correction on the GPU.""" self.update_electric_dev( @@ -549,12 +536,6 @@ class OpenCLPML(PML): self.HPhi1_dev = clarray.to_device(self.queue, self.HPhi1) self.HPhi2_dev = clarray.to_device(self.queue, self.HPhi2) - def set_blocks_per_grid(): - pass - - def get_update_funcs(): - pass - def update_electric(self): """Updates electric field components with the PML correction on the compute device. diff --git a/gprMax/receivers.py b/gprMax/receivers.py index f24c8263..64193636 100644 --- a/gprMax/receivers.py +++ b/gprMax/receivers.py @@ -26,9 +26,7 @@ class Rx: allowableoutputs = ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz", "Ix", "Iy", "Iz"] defaultoutputs = allowableoutputs[:-3] - allowableoutputs_dev = allowableoutputs[:-3] - maxnumoutputs_dev = 0 def __init__(self): self.ID = None @@ -61,9 +59,6 @@ def htod_rx_arrays(G, queue=None): rxcoords[i, 0] = rx.xcoord rxcoords[i, 1] = rx.ycoord rxcoords[i, 2] = rx.zcoord - # Store maximum number of output components - if len(rx.outputs) > Rx.maxnumoutputs_dev: - Rx.maxnumoutputs_dev = len(rx.outputs) # Array to store field components for receivers on compute device - # rows are field components; columns are iterations; pages are receivers diff --git a/gprMax/updates.py b/gprMax/updates.py index 9adece69..c6d103ed 100644 --- a/gprMax/updates.py +++ b/gprMax/updates.py @@ -407,7 +407,9 @@ class CUDAUpdates: knlH = self.source_module(bld, options=config.sim_config.devices["nvcc_opts"]) pml.update_magnetic_dev = knlH.get_function(knl_name) - self._copy_mat_coeffs(knlE, knlH) + # Copy material coefficient arrays to constant memory of GPU - must + # be done for each kernel + self._copy_mat_coeffs(knlE, knlH) def _set_rx_knl(self): """Receivers - initialises arrays on GPU, prepares kernel and gets kernel diff --git a/setup.py b/setup.py index df082040..7c201976 100644 --- a/setup.py +++ b/setup.py @@ -18,6 +18,7 @@ import glob import os +import platform import re import shutil import subprocess @@ -187,35 +188,25 @@ else: ) cpuID = " ".join(cpuID.split()) if "Apple" in cpuID: - gccpath = glob.glob("/opt/homebrew/bin/gcc-[4-9]*") - gccpath += glob.glob("/opt/homebrew/bin/gcc-[10-11]*") - if gccpath: - # Use newest gcc found - os.environ["CC"] = gccpath[-1].split(os.sep)[-1] - rpath = "/opt/homebrew/opt/gcc/lib/gcc/" + gccpath[-1].split(os.sep)[-1][-1] + "/" - else: - raise ( - "Cannot find gcc in /opt/homebrew/bin. gprMax requires gcc " - + "to be installed - easily done through the Homebrew package " - + "manager (http://brew.sh). Note: gcc with OpenMP support " - + "is required." - ) + gccbasepath = "/opt/homebrew/bin/" else: - gccpath = glob.glob("/usr/local/bin/gcc-[4-9]*") - gccpath += glob.glob("/usr/local/bin/gcc-[10-11]*") - if gccpath: - # Use newest gcc found - os.environ["CC"] = gccpath[-1].split(os.sep)[-1] - else: - raise ( - "Cannot find gcc in /usr/local/bin. gprMax requires gcc " - + "to be installed - easily done through the Homebrew package " - + "manager (http://brew.sh). Note: gcc with OpenMP support " - + "is required." - ) + gccbasepath = "/usr/local/bin/" + gccpath = glob.glob(gccbasepath + "gcc-[0-9][0-9]") + if gccpath: + # Use newest gcc found + os.environ["CC"] = gccpath[-1].split(os.sep)[-1] + if "Apple" in cpuID: + rpath = "/opt/homebrew/opt/gcc/lib/gcc/" + gccpath[-1].split(os.sep)[-1][-1] + "/" + else: + raise ( + f"Cannot find gcc in {gccbasepath}. gprMax requires gcc " + + "to be installed - easily done through the Homebrew package " + + "manager (http://brew.sh). Note: gcc with OpenMP support " + + "is required." + ) - # Minimum supported macOS deployment target - MIN_MACOS_VERSION = "10.13" + # Set minimum supported macOS deployment target to installed macOS version + MIN_MACOS_VERSION = platform.mac_ver()[0] try: os.environ["MACOSX_DEPLOYMENT_TARGET"] del os.environ["MACOSX_DEPLOYMENT_TARGET"] diff --git a/testing/analytical_solutions.py b/testing/analytical_solutions.py index 76656dd0..877826d5 100644 --- a/testing/analytical_solutions.py +++ b/testing/analytical_solutions.py @@ -17,8 +17,9 @@ # along with gprMax. If not, see . import numpy as np +from scipy.constants import c +from scipy.constants import epsilon_0 as e0 -import gprMax.config as config from gprMax.waveforms import Waveform @@ -78,35 +79,35 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx): Ex_y = y Ex_z = z - 0.5 * dz Er_x = np.sqrt((Ex_x**2 + Ex_y**2 + Ex_z**2)) - tau_Ex = Er_x / config.sim_config.em_consts["c"] + tau_Ex = Er_x / c # Coordinates of Rx for Ey FDTD component Ey_x = x Ey_y = y + 0.5 * dy Ey_z = z - 0.5 * dz Er_y = np.sqrt((Ey_x**2 + Ey_y**2 + Ey_z**2)) - tau_Ey = Er_y / config.sim_config.em_consts["c"] + tau_Ey = Er_y / c # Coordinates of Rx for Ez FDTD component Ez_x = x Ez_y = y Ez_z = z Er_z = np.sqrt((Ez_x**2 + Ez_y**2 + Ez_z**2)) - tau_Ez = Er_z / config.sim_config.em_consts["c"] + tau_Ez = Er_z / c # Coordinates of Rx for Hx FDTD component Hx_x = x Hx_y = y + 0.5 * dy Hx_z = z Hr_x = np.sqrt((Hx_x**2 + Hx_y**2 + Hx_z**2)) - tau_Hx = Hr_x / config.sim_config.em_consts["c"] + tau_Hx = Hr_x / c # Coordinates of Rx for Hy FDTD component Hy_x = x + 0.5 * dx Hy_y = y Hy_z = z Hr_y = np.sqrt((Hy_x**2 + Hy_y**2 + Hy_z**2)) - tau_Hy = Hr_y / config.sim_config.em_consts["c"] + tau_Hy = Hr_y / c # Initialise fields fields = np.zeros((iterations, 6)) @@ -134,7 +135,7 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx): fdot_Hy = wdot.calculate_value((timestep * dt) - tau_Hy, dt) * dl # Ex - fields[timestep, 0] = ((Ex_x * Ex_z) / (4 * np.pi * config.sim_config.em_consts["e0"] * Er_x**5)) * ( + fields[timestep, 0] = ((Ex_x * Ex_z) / (4 * np.pi * e0 * Er_x**5)) * ( 3 * (fint_Ex + (tau_Ex * f_Ex)) + (tau_Ex**2 * fdot_Ex) ) @@ -145,12 +146,12 @@ def hertzian_dipole_fs(iterations, dt, dxdydz, rx): tmp = 0 fields[timestep, 1] = ( tmp - * ((Ey_x * Ey_z) / (4 * np.pi * config.sim_config.em_consts["e0"] * Er_y**5)) + * ((Ey_x * Ey_z) / (4 * np.pi * e0 * Er_y**5)) * (3 * (fint_Ey + (tau_Ey * f_Ey)) + (tau_Ey**2 * fdot_Ey)) ) # Ez - fields[timestep, 2] = (1 / (4 * np.pi * config.sim_config.em_consts["e0"] * Er_z**5)) * ( + fields[timestep, 2] = (1 / (4 * np.pi * e0 * Er_z**5)) * ( (2 * Ez_z**2 - (Ez_x**2 + Ez_y**2)) * (fint_Ez + (tau_Ez * f_Ez)) - (Ez_x**2 + Ez_y**2) * tau_Ez**2 * fdot_Ez ) diff --git a/testing/diff_output_files.py b/testing/diff_output_files.py new file mode 100644 index 00000000..b2a3da6b --- /dev/null +++ b/testing/diff_output_files.py @@ -0,0 +1,93 @@ +# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom +# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import logging +from pathlib import Path + +import h5py +import numpy as np + +logger = logging.getLogger(__name__) + +def diff_output_files(filename1, filename2): + """Calculates differences between two output files. + + Args: + filename1: string of filename (including path) of output file 1. + filename2: string of filename (including path) of output file 2. + + Returns: + time: numpy array containing time. + datadiffs: numpy array containing power (dB) of differences. + """ + + file1 = h5py.File(Path(filename1), "r") + file2 = h5py.File(Path(filename2), "r") + # Path to receivers in files + path = "rxs/rx1/" + + # Get available field output component names + outputs1 = list(file1[path].keys()) + outputs2 = list(file2[path].keys()) + if outputs1 != outputs2: + logger.exception("Field output components are not the same in each file") + raise ValueError + + # Check that type of float used to store fields matches + floattype1 = file1[path + outputs1[0]].dtype + floattype2 = file2[path + outputs2[0]].dtype + if floattype1 != floattype2: + logger.warning( + f"Type of floating point number in test model ({file1[path + outputs1[0]].dtype}) " + f"does not match type in reference solution ({file2[path + outputs2[0]].dtype})\n" + ) + + # Arrays for storing time + time1 = np.zeros((file1.attrs["Iterations"]), dtype=floattype1) + time1 = np.linspace(0, (file1.attrs["Iterations"] - 1), num=file1.attrs["Iterations"]) + time2 = np.zeros((file2.attrs["Iterations"]), dtype=floattype2) + time2 = np.linspace(0, (file2.attrs["Iterations"] - 1), num=file2.attrs["Iterations"]) + + # Arrays for storing field data + data1 = np.zeros((file1.attrs["Iterations"], len(outputs1)), dtype=floattype1) + data2 = np.zeros((file2.attrs["Iterations"], len(outputs2)), dtype=floattype2) + for ID, name in enumerate(outputs1): + data1[:, ID] = file1[path + str(name)][:] + data2[:, ID] = file2[path + str(name)][:] + if np.any(np.isnan(data1[:, ID])) or np.any(np.isnan(data2[:, ID])): + logger.exception("Data contains NaNs") + raise ValueError + + file1.close() + file2.close() + + # Diffs + datadiffs = np.zeros(data1.shape, dtype=np.float64) + 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 + ) # Replace any division by zero with zero + + # Calculate power (ignore warning from taking a log of any zero values) + with np.errstate(divide="ignore"): + datadiffs[:, i] = 20 * np.log10(datadiffs[:, i]) + # Replace any NaNs or Infs from zero division + datadiffs[:, i][np.invert(np.isfinite(datadiffs[:, i]))] = 0 + + return time1, datadiffs diff --git a/testing/models_pmls/plot_pml_comparison.py b/testing/models_pmls/plot_pml_comparison.py deleted file mode 100644 index 95c5a79f..00000000 --- a/testing/models_pmls/plot_pml_comparison.py +++ /dev/null @@ -1,135 +0,0 @@ -# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom -# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley -# -# This file is part of gprMax. -# -# gprMax is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# gprMax is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with gprMax. If not, see . - -import itertools -import logging -from operator import add -from pathlib import Path - -import h5py -import matplotlib.pyplot as plt -import numpy as np - -logger = logging.getLogger(__name__) - -# Create/setup plot figure -# colors = ['#E60D30', '#5CB7C6', '#A21797', '#A3B347'] # Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php -# colorIDs = ["#62a85b", "#9967c7", "#b3943f", "#6095cd", "#cb5c42", "#c95889"] -colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545", "#ff83ec"] -# colorIDs = ["#ba0044", "#b2d334", "#470055", "#185300", "#ff96b1", "#3e2700", "#0162a9", "#fdb786"] -colors = itertools.cycle(colorIDs) -# for i in range(2): -# next(colors) -lines = itertools.cycle(("--", ":", "-.", "-")) -markers = ["o", "d", "^", "s", "*"] - -parts = Path(__file__).parts -path = "rxs/rx1/" -basename = "pml_3D_pec_plate" -PMLIDs = ["CFS-PML", "HORIPML-1", "HORIPML-2", "MRIPML-1", "MRIPML-2"] -maxerrors = [] -testmodels = ["pml_3D_pec_plate_" + s for s in PMLIDs] - -fig, ax = plt.subplots( - subplot_kw=dict(xlabel="Iterations", ylabel="Error [dB]"), figsize=(20, 10), facecolor="w", edgecolor="w" -) - -for x, model in enumerate(testmodels): - # Open output file and read iterations - fileref = h5py.File(Path(*parts[:-1], basename, basename + "_ref.h5"), "r") - filetest = h5py.File(Path(*parts[:-1], basename, basename + str(x + 1) + ".h5"), "r") - - # Get available field output component names - outputsref = list(fileref[path].keys()) - outputstest = list(filetest[path].keys()) - if outputsref != outputstest: - logger.exception("Field output components do not match reference solution") - raise ValueError - - # Check that type of float used to store fields matches - if filetest[path + outputstest[0]].dtype != fileref[path + outputsref[0]].dtype: - logger.warning( - f"Type of floating point number in test model ({filetest[path + outputstest[0]].dtype}) " - f"does not match type in reference solution ({fileref[path + outputsref[0]].dtype})\n" - ) - floattyperef = fileref[path + outputsref[0]].dtype - floattypetest = filetest[path + outputstest[0]].dtype - # logger.info(f'Data type: {floattypetest}') - - # Arrays for storing time - # timeref = np.zeros((fileref.attrs['Iterations']), dtype=floattyperef) - # timeref = np.linspace(0, (fileref.attrs['Iterations'] - 1) * fileref.attrs['dt'], num=fileref.attrs['Iterations']) / 1e-9 - # timetest = np.zeros((filetest.attrs['Iterations']), dtype=floattypetest) - # timetest = np.linspace(0, (filetest.attrs['Iterations'] - 1) * filetest.attrs['dt'], num=filetest.attrs['Iterations']) / 1e-9 - timeref = np.zeros((fileref.attrs["Iterations"]), dtype=floattyperef) - timeref = np.linspace(0, (fileref.attrs["Iterations"] - 1), num=fileref.attrs["Iterations"]) - timetest = np.zeros((filetest.attrs["Iterations"]), dtype=floattypetest) - timetest = np.linspace(0, (filetest.attrs["Iterations"] - 1), num=filetest.attrs["Iterations"]) - - # Arrays for storing field data - dataref = np.zeros((fileref.attrs["Iterations"], len(outputsref)), dtype=floattyperef) - datatest = np.zeros((filetest.attrs["Iterations"], len(outputstest)), dtype=floattypetest) - for ID, name in enumerate(outputsref): - dataref[:, ID] = fileref[path + str(name)][:] - datatest[:, ID] = filetest[path + str(name)][:] - if np.any(np.isnan(datatest[:, ID])): - logger.exception("Test data contains NaNs") - raise ValueError - - fileref.close() - filetest.close() - - # Diffs - datadiffs = np.zeros(datatest.shape, dtype=np.float64) - for i in range(len(outputstest)): - maxi = np.amax(np.abs(dataref[:, i])) - datadiffs[:, i] = np.divide( - np.abs(datatest[:, i] - dataref[:, i]), maxi, out=np.zeros_like(dataref[:, i]), where=maxi != 0 - ) # Replace any division by zero with zero - - # Calculate power (ignore warning from taking a log of any zero values) - with np.errstate(divide="ignore"): - datadiffs[:, i] = 20 * np.log10(datadiffs[:, i]) - # Replace any NaNs or Infs from zero division - datadiffs[:, i][np.invert(np.isfinite(datadiffs[:, i]))] = 0 - - # Print maximum error value - start = 210 - maxerrors.append(f": {np.amax(datadiffs[start::, 1]):.1f} [dB]") - logger.info(f"{model}: Max. error {maxerrors[x]}") - - # Plot diffs (select column to choose field component, 0-Ex, 1-Ey etc..) - ax.plot(timeref[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)) - ax.set_ylim([-160, -20]) - ax.set_axisbelow(True) - ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed") - -mylegend = list(map(add, PMLIDs, maxerrors)) -legend = ax.legend(mylegend, loc=1, fontsize=14) -frame = legend.get_frame() -frame.set_edgecolor("white") -frame.set_alpha(0) - -plt.show() - -# Save a PDF/PNG of the figure -# fig.savefig(basepath + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1) -# fig.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1) 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 new file mode 100644 index 00000000..52968ef1 --- /dev/null +++ b/testing/models_pmls/pml_3D_pec_plate/plot_pml_comparison.py @@ -0,0 +1,80 @@ +# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom +# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import itertools +import logging +from operator import add +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + +from testing.diff_output_files import diff_output_files + +logger = logging.getLogger(__name__) + +# Create/setup plot figure +# colors = ['#E60D30', '#5CB7C6', '#A21797', '#A3B347'] # Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php +# colorIDs = ["#62a85b", "#9967c7", "#b3943f", "#6095cd", "#cb5c42", "#c95889"] +colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545", "#ff83ec"] +# colorIDs = ["#ba0044", "#b2d334", "#470055", "#185300", "#ff96b1", "#3e2700", "#0162a9", "#fdb786"] +colors = itertools.cycle(colorIDs) +# for i in range(2): +# next(colors) +lines = itertools.cycle(("--", ":", "-.", "-")) +markers = ["o", "d", "^", "s", "*"] + +fn = Path(__file__) +basename = "pml_3D_pec_plate" +PMLIDs = ["CFS-PML", "HORIPML-1", "HORIPML-2", "MRIPML-1", "MRIPML-2"] +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" +) + +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")) + + # Print maximum error value + start = 210 + maxerrors.append(f": {np.amax(datadiffs[start::, 1]):.1f} [dB]") + 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.set_xticks(np.arange(0, 2200, step=100)) + ax.set_xlim([0, 2100]) + ax.set_yticks(np.arange(-160, 0, step=20)) + ax.set_ylim([-160, -20]) + ax.set_axisbelow(True) + ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed") + +mylegend = list(map(add, PMLIDs, maxerrors)) +legend = ax.legend(mylegend, loc=1, fontsize=14) +frame = legend.get_frame() +frame.set_edgecolor("white") +frame.set_alpha(0) + +plt.show() + +# Save a PDF/PNG of the figure +# fig.savefig(basepath + '.pdf', dpi=None, format='pdf', bbox_inches='tight', pad_inches=0.1) +# fig.savefig(savename + '.png', dpi=150, format='png', bbox_inches='tight', pad_inches=0.1) diff --git a/testing/models_pmls/pml_basic/plot_pml_comparison.py b/testing/models_pmls/pml_basic/plot_pml_comparison.py new file mode 100644 index 00000000..86786b07 --- /dev/null +++ b/testing/models_pmls/pml_basic/plot_pml_comparison.py @@ -0,0 +1,82 @@ +# Copyright (C) 2015-2023: The University of Edinburgh, United Kingdom +# Authors: Craig Warren, Antonis Giannopoulos, and John Hartley +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import itertools +import logging +from operator import add +from pathlib import Path + +import matplotlib.gridspec as gridspec +import matplotlib.pyplot as plt +import numpy as np + +from testing.diff_output_files import diff_output_files + +logger = logging.getLogger(__name__) + +# Create/setup plot figure +# Plot colours from http://tools.medialab.sciences-po.fr/iwanthue/index.php +colorIDs = ["#79c72e", "#5774ff", "#ff7c2c", "#4b4e80", "#d7004e", "#007545"] +colors = itertools.cycle(colorIDs) +lines = itertools.cycle(("--", ":", "-.", "-")) +markers = ["o", "d", "^", "s", "*"] + +fn = Path(__file__) +basename = "pml_basic" +PMLIDs = ["off", "x0", "y0", "z0", "xmax", "ymax", "zmax"] +maxerrors = [] + +for x, PMLID in enumerate(PMLIDs): + file1 = fn.parent.joinpath(basename + str(x + 1) + "_CPU.h5") + file2 = fn.parent.joinpath(basename + str(x + 1) + "_GPU.h5") + time, datadiffs = diff_output_files(file1, file2) + + # Print maximum error value + start = 150 + maxerrors.append(f": {np.amax(np.amax(datadiffs[start::, :])):.1f} [dB]") + print(f"{PMLID}: Max. error {maxerrors[x]}") + + # Plot diffs + gs = gridspec.GridSpec(3, 2, hspace=0.3, wspace=0.3) + fig, ax = plt.subplots(figsize=(20, 10), facecolor="w", edgecolor="w") + ax.remove() + fig.suptitle(f"{PMLID}") + + outputs = ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"] + for i, output in enumerate(outputs): + if i < 3: + ax = plt.subplot(gs[i, 0]) + else: + ax = plt.subplot(gs[i - 3, 1]) + ax.plot(time[start::], datadiffs[start::, i], color=next(colors), lw=2) + ax.set_xticks(np.arange(0, 1800, step=200)) + ax.set_xlim([0, 1600]) + ax.set_yticks(np.arange(-400, 80, step=40)) + ax.set_ylim([-400, 40]) + ax.set_axisbelow(True) + ax.grid(color=(0.75, 0.75, 0.75), linestyle="dashed") + ax.set_xlabel("Time [iterations]") + 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 + ".png", dpi=150, format='png', bbox_inches='tight', pad_inches=0.1) + +plt.show() + + diff --git a/testing/models_pmls/pml_basic.py b/testing/models_pmls/pml_basic/pml_basic.py similarity index 100% rename from testing/models_pmls/pml_basic.py rename to testing/models_pmls/pml_basic/pml_basic.py diff --git a/testing/test_models.py b/testing/test_models.py index 196617b9..c4a85348 100644 --- a/testing/test_models.py +++ b/testing/test_models.py @@ -40,14 +40,12 @@ if sys.platform == "linux": python -m testing.test_models """ -# Specify directoty with set of models to test +# Specify directory with set of models to test modelset = "models_basic" # modelset += 'models_advanced' -# modelset += 'models_pmls' basepath = Path(__file__).parents[0] / modelset - # List of available basic test models testmodels = [ "hertzian_dipole_fs_analytical", @@ -64,9 +62,6 @@ testmodels = [ # List of available advanced test models # testmodels = ['antenna_GSSI_1500_fs', 'antenna_MALA_1200_fs'] -# List of available PML models -# testmodels = ['pml_x0', 'pml_y0', 'pml_z0', 'pml_xmax', 'pml_ymax', 'pml_zmax', 'pml_3D_pec_plate'] - # Select a specific model if desired # testmodels = [testmodels[0]] testresults = dict.fromkeys(testmodels) @@ -80,7 +75,7 @@ for i, model in enumerate(testmodels): # Run model file = basepath / model / model - gprMax.run(inputfile=file.with_suffix(".in"), gpu=None) + gprMax.run(inputfile=file.with_suffix(".in"), gpu=None, opencl=None) # Special case for analytical comparison if model == "hertzian_dipole_fs_analytical": @@ -116,6 +111,7 @@ for i, model in enumerate(testmodels): dataref = hertzian_dipole_fs( filetest.attrs["Iterations"], filetest.attrs["dt"], filetest.attrs["dx_dy_dz"], rxposrelative ) + filetest.close() else: # Get output for model and reference files @@ -134,15 +130,15 @@ for i, model in enumerate(testmodels): raise ValueError # Check that type of float used to store fields matches - if filetest[path + outputstest[0]].dtype != fileref[path + outputsref[0]].dtype: - logger.warning( - f"Type of floating point number in test model " - f"({filetest[path + outputstest[0]].dtype}) does not " - f"match type in reference solution ({fileref[path + outputsref[0]].dtype})\n" - ) float_or_doubleref = fileref[path + outputsref[0]].dtype float_or_doubletest = filetest[path + outputstest[0]].dtype - + if float_or_doubleref != float_or_doubletest: + logger.warning( + f"Type of floating point number in test model " + f"({float_or_doubletest}) does not " + f"match type in reference solution ({float_or_doubleref})\n" + ) + # Arrays for storing time timeref = np.zeros((fileref.attrs["Iterations"]), dtype=float_or_doubleref) timeref = ( @@ -166,7 +162,7 @@ for i, model in enumerate(testmodels): raise ValueError fileref.close() - filetest.close() + filetest.close() # Diffs datadiffs = np.zeros(datatest.shape, dtype=np.float64) diff --git a/toolboxes/GPRAntennaModels/GSSI.py b/toolboxes/GPRAntennaModels/GSSI.py index 641cc439..58bf1c9e 100644 --- a/toolboxes/GPRAntennaModels/GSSI.py +++ b/toolboxes/GPRAntennaModels/GSSI.py @@ -420,12 +420,12 @@ def antenna_like_GSSI_1500(x, y, z, resolution=0.001, **kwargs): return scene_objects -def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): +def antenna_like_GSSI_400(x, y, z, resolution=0.002, **kwargs): """Inserts a description of an antenna similar to the GSSI 400MHz antenna. This model represents an update to the previous model of the GSSI 400MHz antenna and was created and optimised by Stadler et al. (2022) in: https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=9686638. - Can be used with 0.5mm, 1mm (default) or 2mm spatial resolution. + Can be used with 2mm spatial resolution. The external dimensions of the antenna are 300x300x178mm. One output point is defined between the arms of the receiver bowtie. The bowties are aligned with the y axis so the output is the y component @@ -451,7 +451,6 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): casesize = (0.3, 0.3, 0.178) # original casethickness = 0.002 shieldthickness = 0.002 - foamsurroundthickness = 0.003 pcbthickness = 0.002 bowtiebase = 0.06 bowtieheight = 0.06 # original 0.056 @@ -488,25 +487,12 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): # Coordinates of source excitation point in antenna tx = x + 0.01 + 0.005 + 0.056, y + casethickness + 0.005 + 0.143, z + skidthickness - if resolution == 0.0005: - dx = 0.0005 - dy = 0.0005 - dz = 0.0005 - tx = x + 0.01 + 0.005 + 0.056, y + casethickness + 0.005 + 0.1435, z + skidthickness - elif resolution == 0.001: - dx = 0.001 - dy = 0.001 - dz = 0.001 - elif resolution == 0.002: - dx = 0.002 - dy = 0.002 - dz = 0.002 - foamsurroundthickness = 0.002 - metalboxheight = 0.088 - tx = x + 0.01 + 0.004 + 0.056, y + casethickness + 0.005 + 0.143 - 0.002, z + skidthickness - 0.002 - else: - logger.exception("This antenna module can only be used with a spatial discretisation of 0.5mm, 1mm, 2mm") - raise ValueError + dx = 0.002 + dy = 0.002 + dz = 0.002 + foamsurroundthickness = 0.002 + metalboxheight = 0.088 + tx = x + 0.01 + 0.004 + 0.056, y + casethickness + 0.005 + 0.143 - 0.002, z + skidthickness - 0.002 # Material definitions absorber = gprMax.Material(er=absorberEr, se=absorbersig, mr=1, sm=0, id="absorber") @@ -565,68 +551,25 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): scene_objects.extend((b1, b2, b3, b4, b5)) # PCB - if resolution == 0.0005: - b6 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.018, y + casethickness + 0.005 + 0.0235, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.034 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - ) - b7 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.178, y + casethickness + 0.005 + 0.0235, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.194 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - ) - scene_objects.extend((b6, b7)) - - elif resolution == 0.001: - b6 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.018, y + casethickness + 0.005 + 0.023, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.034 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - ) - b7 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.178, y + casethickness + 0.005 + 0.023, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.194 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - ) - scene_objects.extend((b6, b7)) - - elif resolution == 0.002: - b6 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), - p2=( - x + 0.01 + 0.005 + 0.033 + bowtiebase, - y + casethickness + 0.006 + 0.202 + patchheight, - z + skidthickness + pcbthickness - 0.002, - ), - material_id="pcb", - ) - b7 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.179, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), - p2=( - x + 0.01 + 0.005 + 0.195 + bowtiebase, - y + casethickness + 0.006 + 0.202 + patchheight, - z + skidthickness + pcbthickness - 0.002, - ), - material_id="pcb", - ) - scene_objects.extend((b6, b7)) + b6 = gprMax.Box( + p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), + p2=( + x + 0.01 + 0.005 + 0.033 + bowtiebase, + y + casethickness + 0.006 + 0.202 + patchheight, + z + skidthickness + pcbthickness - 0.002, + ), + material_id="pcb", + ) + b7 = gprMax.Box( + p1=(x + 0.01 + 0.005 + 0.179, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), + p2=( + x + 0.01 + 0.005 + 0.195 + bowtiebase, + y + casethickness + 0.006 + 0.202 + patchheight, + z + skidthickness + pcbthickness - 0.002, + ), + material_id="pcb", + ) + scene_objects.extend((b6, b7)) elif smooth_dec == "no": # Plastic case @@ -682,366 +625,124 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): scene_objects.extend((b8, b9, b10, b11, b12)) # PCB - if resolution == 0.0005: - b13 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.018, y + casethickness + 0.005 + 0.0235, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.034 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - averaging="n", - ) - b14 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.178, y + casethickness + 0.005 + 0.0235, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.194 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - averaging="n", - ) - scene_objects.extend((b13, b14)) - - elif resolution == 0.001: - b13 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.018, y + casethickness + 0.005 + 0.023, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.034 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - averaging="n", - ) - b14 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.178, y + casethickness + 0.005 + 0.023, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.194 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - averaging="n", - ) - scene_objects.extend((b13, b14)) - - elif resolution == 0.002: - b13 = gprMax.Box( - p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), - p2=( - x + 0.01 + 0.005 + 0.033 + bowtiebase, - y + casethickness + 0.006 + 0.202 + patchheight, - z + skidthickness + pcbthickness - 0.002, - ), - material_id="pcb", - averaging="n", - ) - b14 = gprMax.Box( - 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, - z + skidthickness + pcbthickness, - ), - material_id="pcb", - averaging="n", - ) - scene_objects.extend((b13, b14)) + b13 = gprMax.Box( + p1=(x + 0.01 + 0.005 + 0.017, y + casethickness + 0.005 + 0.021, z + skidthickness - 0.002), + p2=( + x + 0.01 + 0.005 + 0.033 + bowtiebase, + y + casethickness + 0.006 + 0.202 + patchheight, + z + skidthickness + pcbthickness - 0.002, + ), + material_id="pcb", + averaging="n", + ) + b14 = gprMax.Box( + 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, + z + skidthickness + pcbthickness, + ), + material_id="pcb", + averaging="n", + ) + scene_objects.extend((b13, b14)) # PCB components # My own bowties with triangle commands - if resolution == 0.0005: - # "left" side - # extension plates - p1 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.0235, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.026 + bowtiebase, - y + casethickness + 0.005 + 0.0235 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - p2 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.026 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - # triangles - t1 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.0835, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.026 + bowtiebase, y + casethickness + 0.005 + 0.0835, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.026 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.0835 + bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) - t2 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.026 + bowtiebase, y + casethickness + 0.005 + 0.204, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.026 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.204 - bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) + # "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), + p2=( + x + 0.01 + 0.005 + 0.025 + bowtiebase, + y + casethickness + 0.005 + 0.021 + patchheight, + z + skidthickness - 0.002, + ), + material_id="pec", + ) + p2 = gprMax.Plate( + 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, + z + skidthickness - 0.002, + ), + material_id="pec", + ) + # triangles + t1 = gprMax.Triangle( + 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, z + skidthickness - 0.002), + p3=( + x + 0.01 + 0.005 + 0.025 + (bowtiebase / 2), + y + casethickness + 0.005 + 0.081 + bowtieheight, + z + skidthickness - 0.002, + ), + thickness=0, + material_id="pec", + ) + t2 = gprMax.Triangle( + 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, z + skidthickness - 0.002), + p3=( + x + 0.01 + 0.005 + 0.025 + (bowtiebase / 2), + y + casethickness + 0.005 + 0.203 - bowtieheight, + z + skidthickness - 0.002, + ), + thickness=0, + material_id="pec", + ) + # "right" side + p3 = gprMax.Plate( + 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, + z + skidthickness - 0.002, + ), + material_id="pec", + ) + p4 = gprMax.Plate( + 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, + z + skidthickness - 0.002, + ), + material_id="pec", + ) + # triangles + t3 = gprMax.Triangle( + 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, z + skidthickness - 0.002), + p3=( + x + 0.01 + 0.005 + 0.187 + (bowtiebase / 2), + y + casethickness + 0.005 + 0.081 + bowtieheight, + z + skidthickness - 0.002, + ), + thickness=0, + material_id="pec", + ) + t4 = gprMax.Triangle( + 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, z + skidthickness - 0.002), + p3=( + x + 0.01 + 0.005 + 0.187 + (bowtiebase / 2), + y + casethickness + 0.005 + 0.203 - bowtieheight, + z + skidthickness - 0.002, + ), + thickness=0, + material_id="pec", + ) - # "right" side - p3 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.0235, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.186 + bowtiebase, - y + casethickness + 0.005 + 0.0235 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - p4 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.186 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - # triangles - t3 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.0835, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.186 + bowtiebase, y + casethickness + 0.005 + 0.0835, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.186 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.0835 + bowtieheight, - z + skidthickness, - ), - thickness=0, - material_ID="pec", - ) - t4 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.186 + bowtiebase, y + casethickness + 0.005 + 0.204, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.186 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.204 - bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) - - # Edges that represent wire between bowtie halves in 1mm model - e1 = gprMax.Edge(p1=(tx[0] + 0.16, tx[1] - dy, tx[2]), p2=(tx[0] + 0.16, tx[1], tx[2]), material_id="pec") - e2 = gprMax.Edge( - p1=(tx[0] + 0.16, tx[1] + dy, tx[2]), p2=(tx[0] + 0.16, tx[1] + 2 * dy, tx[2]), material_id="pec" - ) - e3 = gprMax.Edge(p1=(tx[0], tx[1] - dy, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec") - e4 = gprMax.Edge(p1=(tx[0], tx[1] + dy, tx[2]), p2=(tx[0], tx[1] + 2 * dy, tx[2]), material_id="pec") - scene_objects.extend((p1, p2, t1, t2, p3, p4, t3, t4, e1, e2, e3, e4)) - - elif resolution == 0.001: - # "left" side - # extension plates - p1 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.023, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.026 + bowtiebase, - y + casethickness + 0.005 + 0.023 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - p2 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.026 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - # triangles - t1 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.083, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.026 + bowtiebase, y + casethickness + 0.005 + 0.083, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.026 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.083 + bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) - t2 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.026, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.026 + bowtiebase, y + casethickness + 0.005 + 0.204, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.026 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.204 - bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) - - # "right" side - p3 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.023, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.186 + bowtiebase, - y + casethickness + 0.005 + 0.023 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - p4 = gprMax.Plate( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=( - x + 0.01 + 0.005 + 0.186 + bowtiebase, - y + casethickness + 0.005 + 0.204 + patchheight, - z + skidthickness, - ), - material_id="pec", - ) - # triangles - t3 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.083, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.186 + bowtiebase, y + casethickness + 0.005 + 0.083, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.186 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.083 + bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) - t4 = gprMax.Triangle( - p1=(x + 0.01 + 0.005 + 0.186, y + casethickness + 0.005 + 0.204, z + skidthickness), - p2=(x + 0.01 + 0.005 + 0.186 + bowtiebase, y + casethickness + 0.005 + 0.204, z + skidthickness), - p3=( - x + 0.01 + 0.005 + 0.186 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.204 - bowtieheight, - z + skidthickness, - ), - thickness=0, - material_id="pec", - ) - - # Edges that represent wire between bowtie halves in 1mm model - e1 = gprMax.Edge(p1=(tx[0] + 0.16, tx[1] - dy, tx[2]), p2=(tx[0] + 0.16, tx[1], tx[2]), material_id="pec") - e2 = gprMax.Edge( - p1=(tx[0] + 0.16, tx[1] + dy, tx[2]), p2=(tx[0] + 0.16, tx[1] + 2 * dy, tx[2]), material_id="pec" - ) - e3 = gprMax.Edge(p1=(tx[0], tx[1] - dy, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec") - e4 = gprMax.Edge(p1=(tx[0], tx[1] + dy, tx[2]), p2=(tx[0], tx[1] + 2 * dy, tx[2]), material_id="pec") - scene_objects.extend((p1, p2, t1, t2, p3, p4, t3, t4, e1, e2, e3, e4)) - - elif resolution == 0.002: - # "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), - p2=( - x + 0.01 + 0.005 + 0.025 + bowtiebase, - y + casethickness + 0.005 + 0.021 + patchheight, - z + skidthickness - 0.002, - ), - material_id="pec", - ) - p2 = gprMax.Plate( - 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, - z + skidthickness - 0.002, - ), - material_id="pec", - ) - # triangles - t1 = gprMax.Triangle( - 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, z + skidthickness - 0.002), - p3=( - x + 0.01 + 0.005 + 0.025 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.081 + bowtieheight, - z + skidthickness - 0.002, - ), - thickness=0, - material_id="pec", - ) - t2 = gprMax.Triangle( - 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, z + skidthickness - 0.002), - p3=( - x + 0.01 + 0.005 + 0.025 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.203 - bowtieheight, - z + skidthickness - 0.002, - ), - thickness=0, - material_id="pec", - ) - # "right" side - p3 = gprMax.Plate( - 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, - z + skidthickness - 0.002, - ), - material_id="pec", - ) - p4 = gprMax.Plate( - 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, - z + skidthickness - 0.002, - ), - material_id="pec", - ) - # triangles - t3 = gprMax.Triangle( - 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, z + skidthickness - 0.002), - p3=( - x + 0.01 + 0.005 + 0.187 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.081 + bowtieheight, - z + skidthickness - 0.002, - ), - thickness=0, - material_id="pec", - ) - t4 = gprMax.Triangle( - 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, z + skidthickness - 0.002), - p3=( - x + 0.01 + 0.005 + 0.187 + (bowtiebase / 2), - y + casethickness + 0.005 + 0.203 - bowtieheight, - z + skidthickness - 0.002, - ), - thickness=0, - material_id="pec", - ) - - # Edges that represent wire between bowtie halves in 2mm model - e1 = gprMax.Edge(p1=(tx[0] + 0.162, tx[1] - dy, tx[2]), p2=(tx[0] + 0.162, tx[1], tx[2]), material_id="pec") - e2 = gprMax.Edge( - p1=(tx[0] + 0.162, tx[1] + dy, tx[2]), p2=(tx[0] + 0.162, tx[1] + 2 * dy, tx[2]), material_id="pec" - ) - e3 = gprMax.Edge(p1=(tx[0], tx[1] - dy, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec") - e4 = gprMax.Edge(p1=(tx[0], tx[1] + dy, tx[2]), p2=(tx[0], tx[1] + 2 * dy, tx[2]), material_id="pec") - scene_objects.extend((p1, p2, t1, t2, p3, p4, t3, t4, e1, e2, e3, e4)) + # Edges that represent wire between bowtie halves in 2mm model + e1 = gprMax.Edge(p1=(tx[0] + 0.162, tx[1] - dy, tx[2]), p2=(tx[0] + 0.162, tx[1], tx[2]), material_id="pec") + e2 = gprMax.Edge( + p1=(tx[0] + 0.162, tx[1] + dy, tx[2]), p2=(tx[0] + 0.162, tx[1] + 2 * dy, tx[2]), material_id="pec" + ) + e3 = gprMax.Edge(p1=(tx[0], tx[1] - dy, tx[2]), p2=(tx[0], tx[1], tx[2]), material_id="pec") + e4 = gprMax.Edge(p1=(tx[0], tx[1] + dy, tx[2]), p2=(tx[0], tx[1] + 2 * dy, tx[2]), material_id="pec") + scene_objects.extend((p1, p2, t1, t2, p3, p4, t3, t4, e1, e2, e3, e4)) # Metallic plate extension b15 = gprMax.Box( @@ -1069,20 +770,18 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): scene_objects.extend((b15, b16)) # Source - # Excitation - Gaussian pulse - w1 = gprMax.Waveform(wave_type="gaussian", amp=1, freq=excitationfreq, id="my_gaussian") - scene_objects.append(w1) - if src_type == "voltage_source": + w1 = gprMax.Waveform(wave_type="gaussian", amp=1, freq=excitationfreq, id="my_gaussian") vs1 = gprMax.VoltageSource( polarisation="y", p1=(tx[0], tx[1], tx[2]), resistance=sourceresistance, waveform_id="my_gaussian" ) - scene_objects.append(vs1) + scene_objects.extend((w1, vs1)) elif src_type == "transmission_line": + w1 = gprMax.Waveform(wave_type="gaussian", amp=1, freq=excitationfreq, id="my_gaussian") tl1 = gprMax.TransmissionLine( polarisation="y", p1=(tx[0], tx[1], tx[2]), resistance=sourceresistance, waveform_id="my_gaussian" ) - scene_objects.append(tl1) + scene_objects.extend((w1, tl1)) else: # Optimised custom pulse exc1 = gprMax.ExcitationFile( @@ -1094,35 +793,19 @@ def antenna_like_GSSI_400(x, y, z, resolution=0.001, **kwargs): scene_objects.extend((exc1, vs1)) # Receiver - # Zero waveform to use with transmission line at receiver output - w2 = gprMax.Waveform(wave_type="gaussian", amp=0, freq=excitationfreq, id="my_zero_wave") - scene_objects.append(w2) - - if resolution == 0.001 or resolution == 0.0005: - if src_type == "transmission_line": - tl2 = gprMax.TransmissionLine( - polarisation="y", - p1=(tx[0] + 0.16, tx[1], tx[2]), - resistance=receiverresistance, - waveform_id="my_zero_wave", - ) - scene_objects.append(tl2) - elif src_type == "voltage_source": - r1 = gprMax.Rx(p1=(tx[0] + 0.16, tx[1], tx[2]), id="rxbowtie", outputs="Ey") - scene_objects.append(r1) - - elif resolution == 0.002: - if src_type == "transmission_line": - tl2 = gprMax.TransmissionLine( - polarisation="y", - p1=(tx[0] + 0.162, tx[1], tx[2]), - resistance=receiverresistance, - waveform_id="my_zero_wave", - ) - scene_objects.append(tl2) - elif src_type == "voltage_source": - r1 = gprMax.Rx(p1=(tx[0] + 0.162, tx[1], tx[2]), id="rxbowtie", outputs="Ey") - scene_objects.append(r1) + if src_type == "transmission_line": + # Zero waveform to use with transmission line at receiver output + w2 = gprMax.Waveform(wave_type="gaussian", amp=0, freq=excitationfreq, id="my_zero_wave") + tl2 = gprMax.TransmissionLine( + polarisation="y", + p1=(tx[0] + 0.162, tx[1], tx[2]), + resistance=receiverresistance, + waveform_id="my_zero_wave", + ) + scene_objects.extend((w2, tl2)) + else: + r1 = gprMax.Rx(p1=(tx[0] + 0.162, tx[1], tx[2]), id="rxbowtie", outputs="Ey") + scene_objects.append(r1) # Geometry views gv1 = gprMax.GeometryView( diff --git a/toolboxes/GPRAntennaModels/README.rst b/toolboxes/GPRAntennaModels/README.rst index a5fcd241..e15b8c0d 100644 --- a/toolboxes/GPRAntennaModels/README.rst +++ b/toolboxes/GPRAntennaModels/README.rst @@ -14,7 +14,7 @@ Manufacturer/Model Dimensions Resolution(s) Author/Contact ======================== ============= ============= ========================================================================================================================================================================================================================= ================ GSSI 1.5GHz (Model 5100) 170x108x45mm 1, 2mm Craig Warren (craig.warren@northumbria.ac.uk), Northumbria University, UK 1,2 MALA 1.2GHz 184x109x46mm 1, 2mm Craig Warren (craig.warren@northumbria.ac.uk), Northumbria University, UK 1 -GSSI 400MHz 300x300x170mm 0.5, 1, 2mm Sam Stadler (Sam.Stadler@liag-hannover.de), `Leibniz Institute for Applied Geophysics `_, Germany 3 +GSSI 400MHz 300x300x170mm 2mm Sam Stadler (sam.stadler@leibniz-liag.de), `Leibniz Institute for Applied Geophysics `_, Germany 3 ======================== ============= ============= ========================================================================================================================================================================================================================= ================ **License**: `Creative Commons Attribution-ShareAlike 4.0 International License `_ @@ -36,7 +36,7 @@ Descriptions of how the models were created can be found in the aforementioned a How to use the package ====================== -The antenna models can be accessed from within a block of Python code in an input file. The models are inserted at location x,y,z. The coordinates are relative to the geometric centre of the antenna in the x-y plane and the bottom of the antenna skid in the z direction. The models must be used with cubic spatial resolutions of either 0.5mm (GSSI 400MHz antenna only), 1mm (default), or 2mm by setting the keyword argument, e.g. ``resolution=0.002``. +The antenna models can be accessed from within a block of Python code in an input file. The models are inserted at location x,y,z. The coordinates are relative to the geometric centre of the antenna in the x-y plane and the bottom of the antenna skid in the z direction. The models must be used with cubic spatial resolutions of either 1mm or 2mm by setting the keyword argument, e.g. ``resolution=0.002``. .. note::