[GSoC 2023] Implementing a discrete plane wave source

[AdittyaPal]

Figure: Waterfall plots of the six auxiliary one-dimensional arrays for the discrete plane wave.

Problem Description

The objective of the project was to implement a discrete plane wave source in gprMax as described in T. Tan and M. Potter, FDTD Discrete Planewave (FDTD-DPW) Formulation for a Perfectly Matched Source in TFSF Simulations, in IEEE Transactions on Antennas and Propagation, vol. 58, no. 8, pp. 2641-2648, Aug. 2010, doi: 10.1109/TAP.2010.2050446. 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 or maybe even multiple targets. Such a source was not formulated in the current version of gprMax used. This pull request is the work done to implement the idea, that has been shown to work in a MATLAB script, and place it in a much larger environment of gprMax, in Python, blending it seamlessly with the other sources available in gprMax.

Deliverables

• added the planeWavesModules.pyx which contains all the cythonized functions necessary for implementing the discrete plane wave in gprMax.
• added class DiscretePlaneWave() to the sources.py which is the blueprint for a discrete plane wave source with the necessary parameters as instance variables.
• added a class DiscretePlaneWave() to cmds_multiuse.py so that objects of the class DiscretePlaneWave() in sources.py can be created from the input provided in the Python API file.
• added to the update_magnetic_sources() method of updates.py to include the plane wave code in a time stepping loop of the FDTD run.
• added a hash command #discrete_plane_wave in hash_cmds_multiuse.py to process the input for a discrete plane wave from a hash input (.in) file.
• added the #discrete_plane_wave hash command to the list of multiplecmds in hash_cmds_file.py.
• added the import statement for the class DiscretePlanewave() so that it can be accessed directly from the Python API in init.py.
• added import statements on the header of the edited files so that the functions can work properly.

Implementation Structure

The plane wave source can be added in the Python API with the creation of an object of the class DiscretePlaneWave() as

plWave = gprMax.DiscretePlaneWave(p1 = (0.010, 0.010, 0.010), p2=(0.040, 0.040, 0.040), psi = 90.0, phi = 63.4, delta_phi = 2.0, theta = 36.7, delta_theta = 1.0, waveform_id = "mypulse")

The parameters p1 and p2 specify the start vertex and the end vertex of the TFSF box, while \psi specifies the polarization, and \phi, \theta specify the orientation of the wave in the TFSF box (in polar coordinates). The waveform_id must match one of the ids of the waveforms initialized. delta_phi, delta_theta are optional parameters to determine how close to the requested phi and theta values the user would like to approximate the rational angles to. There are additional optional parameters, start and stop to control when the place wave source is added and removed from the TFSF box. Lastly, there is a material_id to specify the background medium in the TFSF box where the plane wave propagates. Or alternatively a hash command might be specified

#discrete_plane_wave: 0.010 0.010 0.010 0.040 0.040 0.040 90.0 63.4 36.7 mypulse

with the order of the parameters p1, p2, psi, phi, theta, waveform_id.
Additionally options might be specified in update_magnetic_sources() whether to use the cythonized version of the code or the python version, and whether to precompute the fields values for the sources or compute them on the fly. The workflow roughly follows the scheme:

• A regular three-dimensional FDTD computational domain is set up in the normal fashion using gprMax. The location for the TF/SF interface is determined through the inputted coordinate positions.
• The angle of incidence for the plane-wave is inputted and mapped to integers [ m_x , m_y , m_z ] to get the rational angles by

  \tan \phi = \frac{m_y \Delta x}{m_x \Delta y}
  \tan \theta = \frac{\sqrt{(m_x / \Delta x)^2 + (m_y / \Delta y)^2}}{m_z / \Delta z}

• Six incident plane-wave one-dimensional arrays are initialized to store the 1-D FDTD field grids based on the calculated integers. The source fields for a chunk of cells of these one-dimensional arrays are set up for three of the six fields.
• Normal E field and H field updates are performed in the primary three-dimensional FDTD grid using the existing framework of gprMax.
• Simultaneously, the field updates in the six auxiliary 1-D arrays for the FDTD incident plane-wave using the equations described in the Tan and Potter paper.

The appropriate fields from these auxiliary one-dimensional arrays are placed in on the interface of the total field/scattered field region of the regular three-dimensional FDTD grids in the same iteration of the time stepping loop.

Type of change

• Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected)
• This change requires a documentation update

Checklist:

• I have made this from my own
• I have taken help from some online resources
• My code follows the style guidelines of this project
• I have performed a self-review of my own code
• I have commented my code, particularly in hard-to-understand areas
• I have made corresponding changes to the documentation
• My changes generate no new warnings
• The title of my pull request is a short description of the requested changes.

Future Scope

Currently there is the provision of introducing the plane wave from the direction of the origin only, if the Total field/scattered field box is considered to be present in the first octant (the 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. It might be able to generalize the souring angle in future improvements of the feature.

ATTACH SCREEN-SHOTS

Figure: Progression of the electric fields with time for the discrete plane wave (the fields are plotted in the order E_x, E_y, E_z, constant x section, constant y section, constant z section).

[AdittyaPal]

The changes look sufficient to introduce multiple plane waves in the TF region with minimal leakage in the SF region. I hope it does not break some existing feature of gprMax.