diff --git a/docs/source/_static/css/my_theme.css b/docs/source/_static/css/my_theme.css new file mode 100644 index 00000000..2c85b2c9 --- /dev/null +++ b/docs/source/_static/css/my_theme.css @@ -0,0 +1,27 @@ +@import url("theme.css"); + +/* Main panel */ +.wy-nav-content { + max-width: 80%; +} + + +/* Maths */ +.MathJax_Display { + text-align: left !important; +} + + +/* Citation tables */ +.wy-table-responsive { + margin: 0; + border-top: 1px dotted rgba(235,235,235,1) ! important; +} + +table.citation td { + padding: 5px ! important; +} + +table.citation td.label { + width: 100px ! important; +} diff --git a/docs/source/app_source.rst b/docs/source/app_source.rst new file mode 100644 index 00000000..9b875b97 --- /dev/null +++ b/docs/source/app_source.rst @@ -0,0 +1,202 @@ +*********************** +Overview of source code +*********************** + +This section provides an overview of the source code modules and describes each of the classes and methods used in the gprMax package. The following licensing information applies to all source files unless otherwise stated:: + + Copyright (C) 2015, The University of Edinburgh. + + Authors: Craig Warren and Antonis Giannopoulos + + 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 . + + +constants.py +============ + +Defines constants: + +* Speed of light in vacuum :math:`c=2.9979245 \times 10^8` m/s +* Permittivity of free space :math:`\epsilon_0=8.854187 \times 10^{-12}` F/m +* Permeability of free space :math:`\mu_0=1.256637 \times 10^{-6}` H/m +* Impedance of free space :math:`z_0=376.7303134` Ohms + +Defines data types: + +* Solid and ID arrays use 32-bit integers (0 to 4294967295) +* Rigid arrays use 8-bit integers (the smallest available numpy type to store booleans - true/false) +* Fractal and dispersive coefficient arrays use complex numbers (:code:`complextype`) which are represented as two :code:`floattype` +* Main field arrays use floats (:code:`floattype`) and complex numbers (:code:`complextype`) +* :code:`floattype` and :code:`complextype` are set to use 32-bit floats but can be changed to use 64-bit double precision if required. + +.. automodule:: gprMax.constants + + +exceptions.py +============= + +.. automodule:: gprMax.exceptions + + +fields_output.py +================ + +.. automodule:: gprMax.fields_output + + +fields_update.pyx +================= + +.. automodule:: gprMax.fields_update + + +fractals.py +=========== + +.. automodule:: gprMax.fractals + + +geometry_primitives.pyx +======================= + +.. automodule:: gprMax.geometry_primitives + + +geometry_views.py +================= + +.. automodule:: gprMax.geometry_views + + +gprMax.py +=========== + +.. automodule:: gprMax.gprMax + +grid.py +======= + +.. automodule:: gprMax.grid + + +input_cmds_file.py +================== + +.. automodule:: gprMax.input_cmds_file + + +input_cmds_geometry.py +====================== + +.. automodule:: gprMax.input_cmds_geometry + + +input_cmds_multiuse.py +====================== + +.. automodule:: gprMax.input_cmds_multiuse + + +input_cmds_singleuse.py +======================= + +.. automodule:: gprMax.input_cmds_singleuse + + +materials.py +============ + +.. automodule:: gprMax.materials + + +pml_1order_update.pyx +===================== + +.. automodule:: gprMax.pml_1order_update + + +pml_2order_update.pyx +===================== + +.. automodule:: gprMax.pml_2order_update + + +pml_call_updates.py +=================== + +.. automodule:: gprMax.pml_call_updates + + +pml.py +====== + +.. automodule:: gprMax.pml + + +receivers.py +============ + +.. automodule:: gprMax.receivers + + +snapshots.py +============ + +.. automodule:: gprMax.snapshots + + +sources.py +========== + +.. automodule:: gprMax.sources + + +user_libs.antennas.py +===================== + +This module is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License:: + + Copyright (C) 2015, Craig Warren + + This module is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. + To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/. + + Please use the attribution at http://dx.doi.org/10.1190/1.3548506 + +.. automodule:: user_libs.antennas + + +utilities.py +============ + +.. automodule:: gprMax.utilities + + +waveforms.py +============ + +.. automodule:: gprMax.waveforms + + +yee_cell_build.pyx +================== + +.. automodule:: gprMax.yee_cell_build + + +yee_cell_setget_rigid.pyx +========================= + +.. automodule:: gprMax.yee_cell_setget_rigid + diff --git a/docs/source/app_waveforms.rst b/docs/source/app_waveforms.rst new file mode 100644 index 00000000..babd1aef --- /dev/null +++ b/docs/source/app_waveforms.rst @@ -0,0 +1,154 @@ +.. _waveforms: + +****************** +Built-in waveforms +****************** + +This section provides definitions of the functions that are used to create the built-in waveforms. Example plots are shown using the parameters: amplitude of one, frequency of 1GHz, time window of 6ns, and a time step of 1.926ps. + +gaussian +======== + +A Gaussian waveform. + +.. math:: I = e^{-\zeta(t-\chi)^2} + +where :math:`I` is the current, :math:`\zeta = 2\pi^2f^2`, :math:`\chi=\frac{1}{f}` and :math:`f` is the frequency. + +.. figure:: images/gaussian.pdf + + Example of the ``gaussian`` waveform - time domain and power spectrum. + + +gaussiandot +=========== + +First derivative of a Gaussian waveform. + +.. math:: I = -2 \zeta (t-\chi) e^{-\zeta(t-\chi)^2} + +where :math:`I` is the current, :math:`\zeta = 2\pi^2f^2`, :math:`\chi=\frac{1}{f}` and :math:`f` is the frequency. + +.. figure:: images/gaussiandot.pdf + + Example of the ``gaussiandot`` waveform - time domain and power spectrum. + + +gaussiandotnorm +=============== + +Normalised first derivative of a Gaussian waveform. + +.. math:: I = -2 \sqrt{\frac{e}{2\zeta}} \zeta (t-\chi) e^{-\zeta(t-\chi)^2} + +where :math:`I` is the current, :math:`\zeta = 2\pi^2f^2`, :math:`\chi=\frac{1}{f}` and :math:`f` is the frequency. + +.. figure:: images/gaussiandotnorm.pdf + + Example of the ``gaussiandotnorm`` waveform - time domain and power spectrum. + + +gaussiandotdot +============== + +Second derivative of a Gaussian waveform. + +.. math:: I = 2\zeta \left(2\zeta(t-\chi)^2 - 1 \right) e^{-\zeta(t-\chi)^2} + +where :math:`I` is the current, :math:`\zeta = 2\pi^2f^2`, :math:`\chi=\frac{1}{f}` and :math:`f` is the frequency. + +.. figure:: images/gaussiandotdot.pdf + + Example of the ``gaussiandotdot`` waveform - time domain and power spectrum. + + +gaussiandotdotnorm +================== + +Normalised second derivative of a Gaussian waveform. + +.. math:: I = \left( 2\zeta (t-\chi)^2 - 1 \right) e^{-\zeta(t-\chi)^2} + +where :math:`I` is the current, :math:`\zeta = 2\pi^2f^2`, :math:`\chi=\frac{1}{f}` and :math:`f` is the frequency. + +.. figure:: images/gaussiandotdotnorm.pdf + + Example of the ``gaussiandotdotnorm`` waveform - time domain and power spectrum. + + +gaussiandotdotdot +================= + +Third derivative of a Gaussian waveform. + +.. math:: I = \zeta^2 \left( 3(t-\chi) - 2\zeta (t-\chi)^3 \right) e^{-\zeta(t-\chi)^2} + +where :math:`I` is the current, :math:`\zeta = 2\pi^2f^2`, :math:`\chi=\frac{1}{f}` and :math:`f` is the frequency. + +.. figure:: images/gaussiandotdotdot.pdf + + Example of the ``gaussiandotdotdot`` waveform - time domain and power spectrum. + + +ricker +====== + +A Ricker (or Mexican Hat) waveform which is the negative, normalised second derivative of a Gaussian waveform. + +.. math:: I = - \left( 2\zeta (t-\chi)^2 -1 \right) e^{-\zeta(t-\chi)^2} + +where :math:`I` is the current, :math:`\zeta = 2\pi^2f^2`, :math:`\chi=\frac{1}{f}` and :math:`f` is the frequency. + +.. figure:: images/ricker.pdf + + Example of the ``ricker`` waveform - time domain and power spectrum. + + +sine +==== + +A single cycle of a sine waveform. + +.. math:: I = R\sin(2\pi ft) + +and + +.. math:: + + R = + \begin{cases} + 1 &\text{if $ft\leq1$}, \\ + 0 &\text{if $ft>1$}. + \end{cases} + +:math:`I` is the current, :math:`t` is time and :math:`f` is the frequency. + +.. figure:: images/sine.pdf + + Example of the ``sine`` waveform - time domain and power spectrum. + + +contsine +======== + +A continuous sine waveform. In order to avoid introducing noise into the calculation the amplitude of the waveform is modulated for the first cycle of the sine wave (ramp excitation). + +.. math:: I = R\sin(2\pi ft) + +and + +.. math:: + + R = + \begin{cases} + R_cft &\text{if $R\leq 1$}, \\ + 1 &\text{if $R>1$}. + \end{cases} + +where :math:`I` is the current, :math:`R_c` is set to :math:`0.25`, :math:`t` is time and :math:`f` is the frequency. + +.. figure:: images/contsine.pdf + + Example of the ``contsine`` waveform - time domain and power spectrum. + + diff --git a/docs/source/conf.py b/docs/source/conf.py new file mode 100644 index 00000000..cc0ed043 --- /dev/null +++ b/docs/source/conf.py @@ -0,0 +1,371 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# gprMax documentation build configuration file, created by +# sphinx-quickstart on Fri Jul 17 11:23:46 2015. +# +# This file is execfile()d with the current directory set to its +# containing dir. +# +# Note that not all possible configuration values are present in this +# autogenerated file. +# +# All configuration values have a default; values that are commented out +# serve to show the default. + +import os, re, shlex, sys + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +sys.path.insert(0, os.path.abspath('../../gprMax')) + +# -- General configuration ------------------------------------------------ + +# If your documentation needs a minimal Sphinx version, state it here. +#needs_sphinx = '1.0' + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'sphinx.ext.mathjax', 'sphinx.ext.autodoc', 'sphinx.ext.napoleon' +] + +# Options for autodoc +autodoc_default_flags = ['members'] +autodoc_member_order = 'bysource' +autoclass_content = 'both' + +# Options for napoleon (which parses googlestyle) +napoleon_use_param = False + +# Figure numbering +numfig = True + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# The suffix(es) of source filenames. +# You can specify multiple suffix as a list of string: +# source_suffix = ['.rst', '.md'] +source_suffix = '.rst' + +# The encoding of source files. +#source_encoding = 'utf-8-sig' + +# The master toctree document. +master_doc = 'index' + +# General information about the project. +project = 'gprMax' +copyright = '2015, The University of Edinburgh. Authors: Craig Warren and Antonis Giannopoulos' +author = 'Craig Warren and Antonis Giannopoulos' + +# The version info for the project you're documenting, acts as replacement for +# |version| and |release|, also used in various other places throughout the +# built documents. +# +# The short X.Y version. +# Read version number from gprMax/gprMax.py +#version = re.search('^__version__\s*=\s*\'(.*)\'', +# open('../../gprMax/gprMax/gprMax.py').read(), +# re.M).group(1) +version = '3.0.0b3' +# The full version, including alpha/beta/rc tags. +release = version + '(Bowmore)' + +# The language for content autogenerated by Sphinx. Refer to documentation +# for a list of supported languages. +# +# This is also used if you do content translation via gettext catalogs. +# Usually you set "language" from the command line for these cases. +language = None + +# There are two options for replacing |today|: either, you set today to some +# non-false value, then it is used: +#today = '' +# Else, today_fmt is used as the format for a strftime call. +#today_fmt = '%B %d, %Y' + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +exclude_patterns = ['_build'] + +# The reST default role (used for this markup: `text`) to use for all +# documents. +#default_role = None + +# If true, '()' will be appended to :func: etc. cross-reference text. +#add_function_parentheses = True + +# If true, the current module name will be prepended to all description +# unit titles (such as .. function::). +#add_module_names = True + +# If true, sectionauthor and moduleauthor directives will be shown in the +# output. They are ignored by default. +#show_authors = False + +# The name of the Pygments (syntax highlighting) style to use. +pygments_style = 'sphinx' + +# A list of ignored prefixes for module index sorting. +#modindex_common_prefix = [] + +# If true, keep warnings as "system message" paragraphs in the built documents. +#keep_warnings = False + +# If true, `todo` and `todoList` produce output, else they produce nothing. +todo_include_todos = False + + +# -- Options for HTML output ---------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +html_theme = 'sphinx_rtd_theme' + +# Theme options are theme-specific and customize the look and feel of a theme +# further. For a list of options available for each theme, see the +# documentation. +#html_theme_options = {} + +# Add any paths that contain custom themes here, relative to this directory. +#html_theme_path = [] + +# The name for this set of Sphinx documents. If None, it defaults to +# " v documentation". +html_title = 'gprMax User Guide' + +# A shorter title for the navigation bar. Default is the same as html_title. +#html_short_title = None + +# The name of an image file (relative to this directory) to place at the top +# of the sidebar. +#html_logo = 'images/gprMax_FB_logo.png' + +# The name of an image file (within the static path) to use as favicon of the +# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 +# pixels large. +#html_favicon = None + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] +html_style = 'css/my_theme.css' + +# Add any extra paths that contain custom files (such as robots.txt or +# .htaccess) here, relative to this directory. These files are copied +# directly to the root of the documentation. +#html_extra_path = [] + +# If not '', a 'Last updated on:' timestamp is inserted at every page bottom, +# using the given strftime format. +html_last_updated_fmt = '%b %d, %Y' + +# If true, SmartyPants will be used to convert quotes and dashes to +# typographically correct entities. +#html_use_smartypants = True + +# Custom sidebar templates, maps document names to template names. +#html_sidebars = {} + +# Additional templates that should be rendered to pages, maps page names to +# template names. +#html_additional_pages = {} + +# If false, no module index is generated. +#html_domain_indices = True + +# If false, no index is generated. +#html_use_index = True + +# If true, the index is split into individual pages for each letter. +#html_split_index = False + +# If true, links to the reST sources are added to the pages. +html_show_sourcelink = False + +# If true, "Created using Sphinx" is shown in the HTML footer. Default is True. +#html_show_sphinx = True + +# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True. +#html_show_copyright = True + +# If true, an OpenSearch description file will be output, and all pages will +# contain a tag referring to it. The value of this option must be the +# base URL from which the finished HTML is served. +#html_use_opensearch = '' + +# This is the file name suffix for HTML files (e.g. ".xhtml"). +#html_file_suffix = None + +# Language to be used for generating the HTML full-text search index. +# Sphinx supports the following languages: +# 'da', 'de', 'en', 'es', 'fi', 'fr', 'h', 'it', 'ja' +# 'nl', 'no', 'pt', 'ro', 'r', 'sv', 'tr' +#html_search_language = 'en' + +# A dictionary with options for the search language support, empty by default. +# Now only 'ja' uses this config value +#html_search_options = {'type': 'default'} + +# The name of a javascript file (relative to the configuration directory) that +# implements a search results scorer. If empty, the default will be used. +#html_search_scorer = 'scorer.js' + +# Output file base name for HTML help builder. +htmlhelp_basename = 'gprMaxdoc' + +# -- Options for LaTeX output --------------------------------------------- + +latex_elements = { +# The paper size ('letterpaper' or 'a4paper'). +'papersize': 'a4paper', + +# The font size ('10pt', '11pt' or '12pt'). +#'pointsize': '10pt', + +# Additional stuff for the LaTeX preamble. +#'preamble': '', + +# Latex figure (float) alignment +#'figure_align': 'htbp', +} + +# Grouping the document tree into LaTeX files. List of tuples +# (source start file, target name, title, +# author, documentclass [howto, manual, or own class]). +latex_documents = [ + (master_doc, 'gprMax.tex', 'gprMax User Guide', + 'Craig Warren and Antonis Giannopoulos', 'manual'), +] + +# The name of an image file (relative to this directory) to place at the top of +# the title page. +latex_logo = 'images/gprMax_logo.png' + +# For "manual" documents, if this is true, then toplevel headings are parts, +# not chapters. +#latex_use_parts = False + +# If true, show page references after internal links. +#latex_show_pagerefs = False + +# If true, show URL addresses after external links. +latex_show_urls = 'inline' + +# Documents to append as an appendix to all manuals. +#latex_appendices = [] + +# If false, no module index is generated. +#latex_domain_indices = True + + +# -- Options for manual page output --------------------------------------- + +# One entry per manual page. List of tuples +# (source start file, name, description, authors, manual section). +man_pages = [ + (master_doc, 'gprMax', 'gprMax User Guide', + [author], 1) +] + +# If true, show URL addresses after external links. +#man_show_urls = False + + +# -- Options for Texinfo output ------------------------------------------- + +# Grouping the document tree into Texinfo files. List of tuples +# (source start file, target name, title, author, +# dir menu entry, description, category) +texinfo_documents = [ + (master_doc, 'gprMax', 'gprMax User Guide', + author, 'gprMax', 'One line description of project.', + 'Miscellaneous'), +] + +# Documents to append as an appendix to all manuals. +#texinfo_appendices = [] + +# If false, no module index is generated. +#texinfo_domain_indices = True + +# How to display URL addresses: 'footnote', 'no', or 'inline'. +#texinfo_show_urls = 'footnote' + +# If true, do not generate a @detailmenu in the "Top" node's menu. +#texinfo_no_detailmenu = False + + +# -- Options for Epub output ---------------------------------------------- + +# Bibliographic Dublin Core info. +epub_title = project +epub_author = author +epub_publisher = author +epub_copyright = copyright + +# The basename for the epub file. It defaults to the project name. +#epub_basename = project + +# The HTML theme for the epub output. Since the default themes are not optimized +# for small screen space, using the same theme for HTML and epub output is +# usually not wise. This defaults to 'epub', a theme designed to save visual +# space. +#epub_theme = 'epub' + +# The language of the text. It defaults to the language option +# or 'en' if the language is not set. +#epub_language = '' + +# The scheme of the identifier. Typical schemes are ISBN or URL. +#epub_scheme = '' + +# The unique identifier of the text. This can be a ISBN number +# or the project homepage. +#epub_identifier = '' + +# A unique identification for the text. +#epub_uid = '' + +# A tuple containing the cover image and cover page html template filenames. +#epub_cover = () + +# A sequence of (type, uri, title) tuples for the guide element of content.opf. +#epub_guide = () + +# HTML files that should be inserted before the pages created by sphinx. +# The format is a list of tuples containing the path and title. +#epub_pre_files = [] + +# HTML files shat should be inserted after the pages created by sphinx. +# The format is a list of tuples containing the path and title. +#epub_post_files = [] + +# A list of files that should not be packed into the epub file. +epub_exclude_files = ['search.html'] + +# The depth of the table of contents in toc.ncx. +#epub_tocdepth = 3 + +# Allow duplicate toc entries. +#epub_tocdup = True + +# Choose between 'default' and 'includehidden'. +#epub_tocscope = 'default' + +# Fix unsupported image types using the Pillow. +#epub_fix_images = False + +# Scale large images. +#epub_max_image_width = 0 + +# How to display URL addresses: 'footnote', 'no', or 'inline'. +#epub_show_urls = 'inline' + +# If false, no index is generated. +#epub_use_index = True diff --git a/docs/source/examples_2D.rst b/docs/source/examples_2D.rst new file mode 100644 index 00000000..535d6098 --- /dev/null +++ b/docs/source/examples_2D.rst @@ -0,0 +1,180 @@ +*********** +2D Examples +*********** + +This section provides some general example models in 2D that demonstrate how to use certain features of gprMax. Each example comes with an input file which you can download and run. + +A-scan with a Hertzian dipole source +==================================== + +:download:`cylinder_Ascan_2D.in ` + +This example is the gprMax equivalent of 'Hello World'! It demonstrates how to simulate a single trace (A-scan) from a metal cylinder buried in a dielectric half-space. + +.. literalinclude:: models/cylinder_Ascan_2D.in + :language: none + :linenos: + +The geometry of the scenario is straightforward and is shown in :numref:`cylinder_half_space_geo`. + +.. _cylinder_half_space_geo: + +.. figure:: images/cylinder_half_space_geo.pdf + :width: 600 px + + Geometry of a 2D model of a metal cylinder buried in a dielectric half-space. + +For this initial example a detailed description of what each command in the input file does and why each command was used is given. The following steps explain the steps taken to build the input file: + +Determine the constitutive parameters for the materials +------------------------------------------------------- + +There will be three different materials in the model representing air, the dielectric half-space, and the metal cylinder. Air (free space) already exists as a built-in material in gprMax which can be accessed using the ``free_space`` identifier. The metal cylinder will be modelled as a Perfect Electric Conductor, which again exists as a built-in material in gprMax and can be accessed using the ``pec`` identifier. So the only material which has to be defined is for the dielectric half-space. It is a non-magnetic material, i.e. :math:`\mu_r=1` and :math:`\sigma_*=0` and with a relative permittivity of six, :math:`\epsilon_r=6`, and zero conductivity, :math:`\sigma=0`. The identifier ``half_space`` will be used. + +.. code-block:: none + + #material: 6 0 1 0 half_space + +Determine the source type and excitation frequency +-------------------------------------------------- + +These should generally be known, often based on the GPR system or scenario being modelled. Low frequencies are used where significant penetration depth is important, whereas high frequencies are used where less penetration and better resolution are required. In this case a theoretical Hertzian dipole source fed with a Ricker waveform with a centre frequency of :math:`f_c=1.5~\textrm{GHz}` will be used to simulate the GPR antenna (later examples will demonstrate how to include a model of the actual GPR antenna in the simulation). + +.. code-block:: none + + #waveform: ricker 1 1.5e9 my_ricker + #hertzian_dipole: z 0.100 0.170 0 my_ricker + +The Ricker waveform is created with the ``#waveform`` command, specifying an amplitude of one, centre frequency of 1.5 GHz and picking an arbitrary identifier of ``my_ricker``. The Hertzian dipole source is created using the ``#hertzian_dipole`` command, specifying a z direction polarisation (the survey direction if a B-scan were being created), location on the surface of the slab, and using the Ricker waveform already created. + +Calculate a spatial resolution and domain size +---------------------------------------------- + +In the :ref:`guidance` section it was stated that a good *rule-of-thumb* was that the spatial resolution should be one tenth of the smallest wavelength present in the model. To determine the smallest wavelength, the highest frequency present in the model is required. This is not the centre frequency of the Ricker waveform! By examining the spectrum of the Ricker waveform it is evident much higher frequencies are present, 2-3 times as high as the centre frequency. So the highest frequency present in the model is likely to be around 4 GHz. The wavelength at 4 GHz in the half-space would be: + +.. math:: \lambda = \frac{c}{f \sqrt{\epsilon_r}} = \frac{299792458}{4\times 10^9 \sqrt{6}} = 30.6~\textrm{mm} + +This would give a minimum spatial resolution of 3 mm. However, the diameter of the cylinder is 20 mm so would be resolved to 7 cells. Therefore a better choice would be 2 mm which resolves the diameter of the rebar to 10 cells, but since the model is only 2D (and therefore the size of the domain will not be large), a spatial resolution of 1 mm will be used. + +.. code-block:: none + + #dx_dy_dz: 0.001 0.001 0.001 + +The domain size should be enough to enclose the volume of interest, plus allow 10 cells (if using the default value) for the PML absorbing boundary conditions and approximately another 10 cells of between the PML and any objects of interest. In this case the plan is to take a B-scan of the scenario (in the next example) so the domain should be large enough to do that. Although this is a 2D model one cell must be specified in the infinite direction (in this case the z direction) of the domain. + +.. code-block:: none + + #domain: 0.240 0.190 0.001 + +Since this is a 2D model the PML should be switched off for the two faces of the domain in the infinite direction (in this case the z direction). This is achieved using the command: + +.. code-block:: none + + #pml_cells: 10 10 0 10 10 0 + +Choose a time window +-------------------- + +It is desired to see the reflection from the cylinder, therefore the time window must be long enough to allow the electromagnetic waves to propagate from the source through the half-space to the cylinder and be reflected back to the receiver, i.e. + +.. math:: \frac{0.180}{\frac{c}{\sqrt{6}}} = 1.47~\textrm{ns} + +This is the minimum time required, but to allow for the entire pulse wavelet to be reflected back to the receiver an initial time window of 4 ns will be tested. + +.. code-block:: none + + #time_window: 4e-9 + +gprMax will calculate the time step required for the model using the CFL condition in 3D. However, since this is a 2D model the time step can be relaxed to the CFL condition in 2D. This is achieved using the command: + +.. code-block:: none + + #time_step_limit_type: 2D + +Create the objects +------------------ + +Now physical objects can created for the half-space and the cylinder. First the ``#box`` command will be used to create the half-space and then the ``#cylinder`` command will be given which will overwrite the properties of the half-space with those of the cylinder at the location of the cylinder. + +.. code-block:: none + + #box: 0 0 0 0.240 0.170 0.001 half_space + #cylinder: 0.120 0.080 0 0.120 0.080 0.001 0.010 pec + +Run the model +------------- + +You can now run the model: + +.. code-block:: none + + python -m gprMax cylinder_Ascan_2D.in + +View the results +---------------- + +You should have produced an output file ``cylinder_Ascan_2D.out``. You can view the results (see :ref:`output` section) using the command: + +.. code-block:: none + + python -m tools.plot_hdf5_Ascan cylinder_Ascan_2D.out + +:numref:`cylinder_Ascan_results` shows the time history of the electric and magnetic field components at the receiver location. The :math:`E_z` field component can be converted to voltage which represents the A-scan (trace). The initial part of the signal (<1.5 ns) represents the direct wave from transmitter to receiver. Then comes the reflected wavelet from the metal cylinder. + +.. _cylinder_Ascan_results: + +.. figure:: images/cylinder_Ascan.pdf + + Field outputs from a model of a metal cylinder buried in a dielectric half-space. + + +B-scan with a Hertzian dipole source +==================================== + +:download:`cylinder_Bscan_2D.in ` + +This example using the same geometry as the previous example but this time a B-scan is created. A B-scan is composed of multiple traces (A-scans) recorded as the source and receiver are moved over the target, in this case the metal cylinder. + +.. literalinclude:: models/cylinder_Bscan_2D.in + :language: none + :linenos: + +The differences between this input file and the one from the A-scan are the x coordinates of the source and receiver (lines 11 adnd 12), and the commands needed to move the source and receiver (lines 13 and 14). The source and receiver are offset by 40mm from each other as before. They are now shifted to the starting position for the scan. The ``#src_steps`` command is used to tell gprMax to move every source in the model by specified steps each time the model is run. Similarly, the ``#rx_steps`` command is used to tell gprMax to move every receiver in the model by specified steps each time the model is run. Note, the same functionality could be achieved by using a block of Python code in the input file to move the source and receiver (for further details see the :ref:`Python section `). + +To run the model for a B-scan is slightly different than for a single A-scan. You must pass an optional argument to gprMax to specify the number of times the model should be run, which in this case is the number of A-scans (traces) that will comprise the B-scan. For a B-scan over a distance of 160mm with a step of 2mm that is 80 A-scans. + +.. code-block:: none + + python -m gprMax cylinder_Bscan_2D.in -n 80 + + +Results +------- + +You should have produced 80 output files, one for each A-scan, with names ``cylinder_Bscan_2D1.out``, ``cylinder_Bscan_2D2.out`` etc... These can be combined into a single file using the command: + +.. code-block:: none + + python -m tools.outputfiles_mergehdf5 cylinder_Bscan_2D 80 + +You should see a combined output file ``cylinder_Bscan_2D_all.out``. The tool will ask you if you want to delete the original single A-scan output files or keep them. + +You can now view an image of the B-scan using the command: + +.. code-block:: none + + python -m tools.plot_hdf5_Bscan cylinder_Bscan_2D_all.out Ez + +:numref:`cylinder_Bscan_results` shows the B-scan (image of the Ez field). As expected a hyperbolic response is present from the metal cylinder. + +.. _cylinder_Bscan_results: + +.. figure:: images/cylinder_Bscan_results.pdf + + B-scan of model of a metal cylinder buried in a dielectric half-space. + + + + + + diff --git a/docs/source/examples_3D.rst b/docs/source/examples_3D.rst new file mode 100644 index 00000000..4ce4f4c3 --- /dev/null +++ b/docs/source/examples_3D.rst @@ -0,0 +1,110 @@ +*********** +3D Examples +*********** + +This section provides some general example models in 3D that demonstrate how to use certain features of gprMax. Each example comes with an input file which you can download and run. + +Using an antenna model +====================== + +:download:`antenna_MALA_1200_fs.in ` + +This example demonstrates how to use one of the built-in antenna models in a simulation. Using a model of an antenna rather than a simple source, such as a Hertzian dipole, can improve the accuracy of the results of a simulation for many situations. It is especially important when the target is in the near-field of the antenna and there are complex interactions between the antenna and the environment. The simulation uses the model of an antenna similar to a MALA 1.2GHz antenna. + +.. literalinclude:: models/antenna_MALA_1200_fs.in + :language: none + :linenos: + +.. figure:: images/antenna_MALA_1200.png + :width: 600 px + + FDTD geometry mesh showing an antenna model similar to a MALA 1.2GHz antenna (skid removed for illustrative purposes). + +The antenna model is loaded from a Python module and inserted into the input file just like another geometry command. The arguments for the ``antenna_like_MALA_1200`` function specify its (x, y, z) location as 0.132m, 0.095m, 0.100m using a 1mm spatial resolution. In this example the antenna is the only object in the model, i.e. the antenna is in free space. When the simulation is run two geometry files for the antenna are produced along with an output file which contains a single receiver (the antenna output). The antenna bowties are aligned with the y axis so the output is the y component of the electric field. More information can be found in the :ref:`Python section `. + +Results +------- + +:numref:`antenna_MALA_1200_fs_results` shows the time history of the electric and magnetic field components from the receiver bowtie of the antenna model. The antenna bowties are aligned with the y axis so the output will be the Ey component of the electric field. + +.. _antenna_MALA_1200_fs_results: + +.. figure:: images/antenna_MALA_1200_fs_results.pdf + + Field outputs from the receiver bowtie of a model of an antenna similar to a MALA 1.2GHz antenna. + +B-scan with an antenna model +============================ + +:download:`GSSI_1500_cylinder_Bscan.in ` + +This example demonstrates how to create a B-scan with an antenna model. The scenario is purposely simple to illustrate the method. A metal cylinder of diameter 20mm is buried in a dielectric half-space which has a relative permittivity of six. The simulation uses the model of an antenna similar to a GSSI 1.5GHz antenna. + +.. literalinclude:: models/GSSI_1500_cylinder_Bscan.in + :language: none + :linenos: + +.. figure:: images/GSSI_1500_cylinder.png + :width: 600 px + + FDTD geometry mesh showing a metal cylinder buried in a half-space and an antenna model similar to a GSSI 1.5GHz antenna. + +The antenna must be moved to a new position for every single A-scan (trace) in the B-scan. In this example the B-scan distance will be 270mm with a trace every 5mm, so 54 model runs will be required. + +.. code-block:: none + + python -m gprMax GSSI_1500_cylinder_Bscan.in -n 54 + +The total number of runs for a model as well as the number of the current run of the model are stored and can be accessed in Python as ``number_model_runs`` and ``current_model_run``. The ``current_model_run`` can be used to move the position of the antenna for every run of the model as shown in Line 13. The antenna will be moved 5mm in the x direction for every new run of the model. + +Results +------- + +:numref:`GSSI_1500_cylinder_Bscan_results` shows the B-scan (image of the Ey field). As expected a hyperbolic response is present from the metal cylinder. + +.. _GSSI_1500_cylinder_Bscan_results: + +.. figure:: images/GSSI_1500_cylinder_Bscan_results.pdf + + B-scan of model of a metal cylinder buried in a dielectric half-space with a model of an antenna similar to a GSSI 1.5GHz antenna. + +Building a heterogeneous soil +============================= + +:download:`heterogeneous_soil.in ` + +This example demonstrates how to build a more realistic soil model using a stochastic distribution of dielectric properties. A mixing model for soils proposed by Peplinski (http://dx.doi.org/10.1109/36.387598) is used to define a series of dispersive material properties for the soil. + +.. literalinclude:: models/heterogeneous_soil.in + :language: none + :linenos: + +.. figure:: images/heterogeneous_soil.png + :width: 600 px + + FDTD geometry mesh showing a heterogeneous soil model with a rough surface. + +Line 10 defines a series of dispersive materials to represent a soil with sand fraction 0.5, clay fraction 0.5, bulk density :math:`2~g/cm^3`, sand particle density of :math:`2.66~g/cm^3`, and a volumetric water fraction range of 0.001 - 0.25. The volumetric water fraction is given as a range which is what defines a series of dispersive materials. + +These materials can then be distributed stochastically over a volume using the ``#fractal_box`` command. Line 11 defines a volume, a fractal dimension, a number of materials, and a mixing model to use. The fractal dimension, 1.5, controls how the materials are stochastically distributed. The fractal weightings, 1, 1, 1, weight the fractal in the x, y, and z directions. The number of materials, 50, specifies how many dispersive materials to create using the mixing model (``my_soil``). + +Adding rough surfaces +--------------------- + +A rough surface can be added to any side of ``#fractal_box`` using, + +.. code-block:: none + + #add_surface_roughness: 0 0 0.070 0.15 0.15 0.070 1.5 1 1 0.065 0.080 my_soil_box + +which defines one of the surfaces of the ``#fractal_box``, a fractal dimension, and minimum and maximum values for the height of the roughness (relative to the original ``#fractal_box`` volume). In this example the roughness will be stochastically distributed with troughs up to 5mm deep, and peaks up to 10mm high. + +More information, including adding surface water and vegetation, can be found in the :ref:`section on using the fractal box command `. + + + + + + + + diff --git a/docs/source/faqs.rst b/docs/source/faqs.rst new file mode 100644 index 00000000..8259835b --- /dev/null +++ b/docs/source/faqs.rst @@ -0,0 +1,5 @@ +**** +FAQs +**** + +This section will provide answers to frequently asked questions about gprMax and its uses. \ No newline at end of file diff --git a/docs/source/features.rst b/docs/source/features.rst new file mode 100644 index 00000000..f815a502 --- /dev/null +++ b/docs/source/features.rst @@ -0,0 +1,163 @@ +.. _capabilities: + +***************** +Software Features +***************** + +This section begins with an overview, primarily for previous users, about what is new, what has changed, and what has been retired in this version of gprMax. It is then followed by more general descriptions some of the key features of gprMax that are useful for GPR modelling as well as more general electromagnetic simulations. + +What's new/changed? +=================== + +A brief summary of what each input command does is given. Please refer to the :ref:`commands` section for a detailed description of the syntax of each command. + +The code has been completely re-written in Python/Cython. In the process a lot of changes have been made to improve efficiency, speed, and usability. Many new features have been implemented, which have been focussed on the following areas: + +* Scripting in the input file. +* Built-in library of antenna models +* Anisotropic material modelling. +* Dispersive material modelling using multiple pole Debye, Lorenz or Drude formulations. +* Building heterogeneous objects using fractal distributions. +* Building objects with rough surfaces. +* Modelling soils with realistic dielectric and geometric properties. +* Improved PML (RIPML) performance. + +New commands +------------ + +* ``#python`` and ``#end_python`` are used to define blocks of the input file where Python code will be executed. This allows the user to use scripting directly in the input file. +* ``#material`` replaces ``#medium`` with a new syntax. +* ``#add_dispersion_debye`` is used to add Debye dispersive properties to a ``#material``. +* ``#add_dispersion_lorenz`` is used to add Lorenz dispersive properties to a ``#material``. +* ``#add_dispersion_drude`` is used to add Drude dispersive properties to a ``#material``. +* ``#soil_peplinski`` is a soil mixing model that can be used with ``#fractal_box`` to generate soil(s) with more realistic dielectric and geometric properties. +* ``#cylindrical_sector`` (like a slice of pie shape) is a new object building command. +* ``#geometry_view`` replaces ``#geometry_file`` or ``#geometry_vtk`` and is used to create views of the geometry of the model in open source Visualization ToolKit (VTK) (http://www.vtk.org) format which can be viewed in many free readers, such as Paraview (http://www.paraview.org). +* ``#fractal_box`` is used to create a volume with a fractal distribution of properties. +* ``#add_surface_roughness`` is used to add a rough surface to a ``#fractal_box``. +* ``#add_surface_water`` is used to add surface water to a ``#fractal_box`` that has a rough surface. +* ``#add_grass`` is used to add grass to a ``#fractal_box``. +* ``#waveform`` is used to specify a waveform shape, and works in conjunction with the changed source commands. +* ``#magnetic_dipole`` is used to introduce a magnetic dipole source, i.e. current on a small loop. +* ``#pml_cfs`` is an advanced command for adjusting the CFS parameters used in the new PML. + + +Changed commands +---------------- + +* All object building commands support anisotropy via additional material identifiers, e.g. ``#box: 0.0 0.0 0.0 0.1 0.1 0.1 matX matY matZ``. +* Dielectric smoothing can be turned on (the default setting) or off for any volumetric object building command by specifying a character ``y`` (on) or ``n`` (off) after the material identifier, e.g. ``#sphere: 0.5 0.5 0.5 0.25 sand n``. +* ``#triangle`` can now create both triangular patches (2D) and triangular prisms (3D) via a thickness parameter. +* ``#pml_cells`` replaces ``#pml_layers`` and can now be used to control the number of cells of PML on the six faces of the model domain. The number of cells can be set to zero on any of the faces to turn that PML off if desired. The default behaviour (if this command is not specified) is to use 10 cells. +* ``#hertzian_dipole`` and ``#voltage_source`` now specify polarisation, location, any additional parameters, and an identifier to link to a ``#waveform`` command. +* ``#snapshot`` no longer requires a type, as only VTK snapshot files are now produced. +* ``#num_of_procs`` is now called ``#num_omp_threads``. +* ``#tx_steps`` is now called ``#src_steps``. + + +Retired commands +---------------- + +* ``#analysis`` and ``#end_analysis`` are no longer required as: sources and receivers can be specified anywhere is the input file; the output file automatically has the same name as the input file with a ``.out`` extension; the format of the output file is now HDF5 (https://www.hdfgroup.org/HDF5/); gprMax can be run with the syntax ``gprMax my_input_file -n number_of_runs``, where ``number_of_runs`` can be used to rerun the model for creating scans and/or moving geometry between runs. +* ``#tx`` is no longer required as the polarisation and position of a source is now specified in the source command, e.g. ``#hertzian_dipole: y 0.05 0.05 0.05 myPulse``. +* ``#cylinder_new`` has become ``#cylinder``. +* ``#cylindrical_segment`` was under-used and its effect can be created by cutting a ``#cylinder`` with a ``#box``. +* ``#bowtie`` can be created using the new behaviour of ``#triangle``. +* ``#number_of_media`` is not required with the new Python code. +* ``#nips_number`` is not required with the new Python code. +* ``#media_file`` was under-used. +* ``#geometry_file`` has been replaced with the ``#geometry_view`` command. +* ``#medium`` has been replaced with the ``#material`` command. +* ``#abc_type``, ``#abc_order``, ``#abc_stability_factors``, ``#abc_optimization_angles``, and ``#abc_mixing_parameters`` are not required as the Higdon ABCs have been removed. +* ``#huygens_surface`` will become part of the new ``#plane_wave`` command which is yet to be implemented. + + +Commands yet to be implemented +------------------------------ + +There are commands from previous versions of gprMax that are planned for this version, but are yet to be implemented. These will be introduced in a future update. They are: ``#thin_wire``, ``#transmission_line``, and ``#plane_wave``. + + +Migrating old input files +------------------------- + +gprMax includes a Python module (in the ``tools`` package) to help you migrate old input files, written for previous versions of gprMax, to the syntax of the new commands. The module will do its best to convert the old file and write a new one, however, you should still carefully check the new file to make sure it is what you intended! Usage (from the top-level gprMax directory) is: ``python -m tools.inputfileold2new my_old_inputfile.in``. + + +Key features +============ + +Python scriptable input files +----------------------------- + +The input file has now been made scriptable by permitting blocks of Python code to be specified between ``#python`` and ``#end_python`` commands. The code is executed when the input file is read by gprMax. You don't need any external tools, such as MATLAB, to generate larger, more complex input files for building intricate models. Python scripting means that gprMax now includes :ref:`libraries of more complex objects, such as antennas `, that can be easily inserted into a model. You can also access a number of built-in constants from your Python code. For further details see the :ref:`Python section `. + +Dispersive media +---------------- + +gprMax has always included the ability to represent dispersive materials using a single-pole Debye model. Many materials can be adequately represented using this approach for the typical frequency ranges associated with GPR. However, multi-pole Debye, Drude and Lorenz functions are often used to simulate the electric susceptibility of materials such as: water [PIE2009]_, human tissue [IRE2013]_, cold plasma [LI2013]_, gold [VIA2005]_, and soils [BER1998]_, [GIAK2012]_, [TEI1998]_. Electric susceptibility relates the polarization density to the electric field, and includes both the real and imaginary parts of the complex electric permittivity variation. In the new version of gprMax a recursive convolution based method is used to express dispersive properties as apparent current density sources [GIA2014]_. A major advantage of this implementation is that it creates an inclusive susceptibility function that holds, as special cases, Debye, Drude and Lorenz materials. For further details see the :ref:`material commands section `. + +Realistic soils, heterogeneous objects and rough surfaces +--------------------------------------------------------- + +The inclusion of improved models of soils is important for many GPR simulations. gprMax can now be used to create soils with more realistic dielectric and geometrical properties. A semi-empirical model, initially suggested by [DOB1985]_, is used to describe the dielectric properties of the soil. The model relates relative permittivity of the soil to bulk density, sand particle density, sand fraction, clay fraction and water volumetric fraction. Using this approach, a more realistic soil with a stochastic distribution of the aforementioned parameters can be modelled. The real and imaginary parts of this semi-empirical model can be approximated using a multi-pole Debye function plus a conductive term. This can now be achieved in gprMax using the new dispersive material functionality. For further details see the :ref:`material commands section `. + +Fractals are scale invariant functions which can express the topography of the earth for a wide range of scales with sufficient detail [TUR1987]_. For this reason fractals have been chosen to represent the topography of soils. Fractals can be generated by the convolution of Gaussian noise with an inverse Fourier transform of :math:`\frac{1}{kb}`, where :math:`k` is the wavenumber and :math:`b` is a constant related to the fractal dimension [TUR1997]_. gprMax can now generate heterogeneous volumes (boxes) with realistic soil properties that can have rough surfaces applied. For further details see the :ref:`fractal object building commands section `. + +Fractal correlated noise [TUR1997]_ is used to describe the stochastic distribution of the properties of soils. This approach has been chosen because it has been shown that soil-related environmental properties frequently obey fractal laws [BUR1981]_, [HILL1998]_. For further details see the :ref:`material commands section ` and the :ref:`fractal object building commands section `. + +.. _antennas: + +Library of antenna models +------------------------- + +gprMax now includes Python modules with pre-defined models of antennas that behave similarly to commercial antennas [WAR2011]_. Currently models of antennas similar to Geophysical Survey Systems, Inc. (GSSI) (http://www.geophysical.com) 1.5 GHz (Model 5100) antenna, and MALA Geoscience (http://www.malags.com/) 1.2 GHz antenna are included. By taking advantage of Python scripting in input files, using such complex structures in a model is straightforward without having to be built step-by-step by the user. For further details see the :ref:`Python section `. + +Anisotropy +---------- + +It is possible to specify objects that have diagonal anisotropy which allows materials such as wood and fibre-reinforced composites, often imaged with GPR, to be more accurately modelled. Standard isotropic objects specify one material identifier that defines the same properties in x, y, and z directions. However, every volumetric object building command can also be specified with three material identifiers, which allows properties for the x, y, and z directions to be separately defined. + +Dielectric smoothing +-------------------- + +At the boundaries between different materials in the model there is the question of which material properties to use? + +* Should the last object to be defined at that location dictate the properties? +* Should an average set of properties of the materials of the objects that share that location be used? + +This latter option is often referred to as dielectric smoothing. To address this question gprMax includes an option to turn dielectric smoothing on or off for volumetric object building commands. The default behaviour (if no option is specified) is for dielectric smoothing to be on. The option can be specified with a single character ``y`` (on) or ``n`` (off) given after the material identifier in each object command. + +Perfectly Matched Layer (PML) boundary conditions +------------------------------------------------- + +With increased research into quantitative information from GPR, it has become necessary for models to be able to have more efficient and better-performing Perfectly Matched Layer (PML) absorbing boundary conditions. Since 2005 gprMax has featured PML absorbing boundary conditions based on the uniaxial PML (UPML) [GED1998]_ formulation. A PML based on a recursive integration approach to the complex frequency shifted (CFS) PML [GIA2012]_ has been adopted in the new version of gprMax. A general formulation of this RIPML, which can be used to develop any order of PML, has been used to implement first and second order CFS stretching functions. One of the attractions of the RIPML is that it is easily applied as a correction to the field quantities after the complete FDTD grid has been updated using the standard FDTD update equations. gprMax now offers the ability (for advanced users) to customise the parameters of the PML which allows its performance to be better optimised for specific applications. Additionally, since the RIPML is media agnostic it can be used without change to problems involving dispersive and anisotropic materials. For further details see the :ref:`PML commands section `. + +Open source, robust, file formats +--------------------------------- + +Alongside improvements to the input file there is a new output file format – HDF5 (http://www.hdfgroup.org/HDF5/) – to manage the larger and more complex data sets that are being generated. HDF5 is a robust, portable and extensible format with a number of free readers available. For further details see the :ref:`output file section `. + +In addition, the Visualization Toolkit (VTK) (http://www.vtk.org) is being used for improved handling and viewing of the detailed 3D FDTD geometry meshes. The VTK is an open-source system for 3D computer graphics, image processing and visualisation. It also has a number of free readers available including Paraview (http://www.paraview.org). For further details see the :ref:`geometry view command `. + +**References** + +.. [PIE2009] Pieraccini, M., Bicci, A., Mecatti, D., Macaluso, G., & Atzeni, C. (2009). Propagation of large bandwidth microwave signals in water. Antennas and Propagation, IEEE Transactions on, 57(11), 3612-3618. (http://dx.doi.org/10.1109/tap.2009.2025674) +.. [IRE2013] Ireland, D., & Abbosh, A. (2013). Modeling human head at microwave frequencies using optimized Debye models and FDTD method. Antennas and Propagation, IEEE Transactions on, 61(4), 2352-2355. (http://dx.doi.org/10.1109/tap.2013.2242037) +.. [LI2013] Li, J., Guo, L. X., Jiao, Y. C., & Wang, R. (2013). Composite scattering of a plasma-coated target above dispersive sea surface by the ADE-FDTD method. Geoscience and Remote Sensing Letters, IEEE, 10(1), 4-8. (http://dx.doi.org/10.1109/lgrs.2012.2189751) +.. [VIA2005] Vial, A., Grimault, A. S., Macías, D., Barchiesi, D., & de La Chapelle, M. L. (2005). Improved analytical fit of gold dispersion: Application to the modeling of extinction spectra with a finite-difference time-domain method. Physical Review B, 71(8), 085416. (http://dx.doi.org/10.1103/physrevb.71.085416) +.. [BER1998] Bergmann, T., Robertsson, J. O., & Holliger, K. (1998). Finite-difference modeling of electromagnetic wave propagation in dispersive and attenuating media. Geophysics, 63(3), 856-867. (http://dx.doi.org/10.1190/1.1444396) +.. [GIAK2012] Giannakis, I., Giannopoulos, A., & Davidson, N. (2012). Incorporating dispersive electrical properties in FDTD GPR models using a general Cole-Cole dispersion function. In 2012 14th International Conference on Ground Penetrating Radar (GPR). (http://dx.doi.org/10.1109/icgpr.2012.6254866) +.. [TEI1998] Teixeira, F. L., Chew, W. C., Straka, M., Oristaglio, M. L., & Wang, T. (1998). Finite-difference time-domain simulation of ground penetrating radar on dispersive, inhomogeneous, and conductive soils. Geoscience and Remote Sensing, IEEE Transactions on, 36(6), 1928-1937. (http://dx.doi.org/10.1109/36.729364) +.. [GIA2014] Giannakis, I., & Giannopoulos, A. (2014). A Novel Piecewise Linear Recursive Convolution Approach for Dispersive Media Using the Finite-Difference Time-Domain Method. Antennas and Propagation, IEEE Transactions on, 62(5), 2669-2678. (http://dx.doi.org/10.1109/tap.2014.2308549) +.. [DOB1985] Dobson, M. C., Ulaby, F. T., Hallikainen, M. T., & El-Rayes, M. (1985). Microwave dielectric behavior of wet soil-Part II: Dielectric mixing models. Geoscience and Remote Sensing, IEEE Transactions on, (1), 35-46. (http://dx.doi.org/10.1109/tgrs.1985.289498) +.. [TUR1987] Turcotte, D. L. (1987). A fractal interpretation of topography and geoid spectra on the Earth, Moon, Venus, and Mars. Journal of Geophysical Research: Solid Earth (1978–2012), 92(B4), E597-E601. (http://dx.doi.org/10.1029/jb092ib04p0e597) +.. [TUR1997] Turcotte, D. L. (1997). Fractals and chaos in geology and geophysics. Cambridge university press. (http://dx.doi.org/10.1017/cbo9781139174695) +.. [BUR1981] Burrough, P. A. (1981). Fractal dimensions of landscapes and other environmental data. Nature, 294(5838), 240-242. (http://dx.doi.org/10.1038/294240a0) +.. [HILL1998] Hillel, D. (1998). Environmental soil physics: Fundamentals, applications, and environmental considerations. Academic press. (http://dx.doi.org/10.1016/b978-012348525-0/50030-6) +.. [WAR2011] Warren, C., & Giannopoulos, A. (2011). Creating finite-difference time-domain models of commercial ground-penetrating radar antennas using Taguchi’s optimization method. Geophysics, 76(2), G37-G47. (http://dx.doi.org/10.1190/1.3548506) +.. [GED1998] Gedney, S. D. (1998). The perfectly matched layer absorbing medium. Advances in Computational Electrodynamics: The Finite-Difference Time-Domain Method, 263-344. +.. [GIA2012] Giannopoulos, A. (2012). Unsplit implementation of higher order PMLs. Antennas and Propagation, IEEE Transactions on, 60(3), 1479-1485. (http://dx.doi.org/10.1109/tap.2011.2180344) + + + diff --git a/docs/source/geometry_snapshots.rst b/docs/source/geometry_snapshots.rst new file mode 100644 index 00000000..f156b740 --- /dev/null +++ b/docs/source/geometry_snapshots.rst @@ -0,0 +1,53 @@ +*************************** +Geometry and Snapshot files +*************************** + +The geometry and snapshot files use the open source Visualization ToolKit (VTK) (http://www.vtk.org) format which can be viewed in many free readers, such as Paraview (http://www.paraview.org). Paraview is an open-source, multi-platform data analysis and visualization application. It is available for Linux, Mac OS X, and Windows. + +Geometry files +============== + +The ``#geometry_view:`` command produces either ImageData (.vti) for a per-cell geometry view, or PolygonalData (.vtp) for a per-cell-edge geometry view. The following are steps to get started with viewing geometry files in Paraview: + +.. _pv_toolbar: + +.. figure:: images/paraview_toolbar.png + + Paraview toolbar showing ``gprMax_materials`` macro button. + +#. **Open the file** either from the File menu or toolbar. +#. Click the **Apply** button in the Properties panel. You should see an outline of the volume of the geometry view. +#. Install the ``gprMax_materials.py`` Python script, that comes with the gprMax source code (in the ``tools/Paraview macros`` directory), as a macro in Paraview. This script makes it quick and easy to view the different materials in a geometry file. To add the script as a macro in Paraview choose the file from the Macros->Add new macro menu. It will then appear as a shortcut button in the toolbar as shown in :numref:`pv_toolbar`. You only need to do this once, the macro will be kept in Paraview for future use. +#. Click the ``gprMax_materials`` shortcut button. All the materials in the model should appear in the Pipeline Browser as Threshold items as shown in :numref:`pv_pipeline`. + +.. _pv_pipeline: + +.. figure:: images/paraview_pipeline.png + :width: 350 px + + Paraview Pipeline Browser showing list of materials in an example model. + +.. tip:: + * You can turn on and off the visibility of materials using the eye icon in the Pipeline Browser. You can select multiple materials using the Shift key, and by shift-clicking the eye icon, turn the visibility of multiple materials on and off. + + * You can set the Color and Opacity of materials from the Properties panel. + + +Snapshot files +============== + +The ``#snapshot:`` command produces an ImageData (.vti) snapshot file for each time instance requested. The following are steps to get started with viewing snapshot files in Paraview: + +#. **Open the file** either from the File menu or toolbar. Paraview should recognise the time series based on the file name and load in all the files. +#. Click the **Apply** button in the Properties panel. You should see an outline of the snapshot volume. +#. Use the **Coloring** drop down menu to select either **E-field** or **H-field**, and the further drop down menu to select either **Magnitude**, **x**, **y** or **z** component. +#. From the **Representation** drop down menu select **Surface**. +#. You can step through or play as an animation the time steps using the **time controls** in the toolbar. + +.. tip:: + + * Turn on the Animation View (View->Animation View menu) to control the speed and start/stop points of the animation. + + * Use the Color Map Editor to adjust the Color Scaling. + + diff --git a/docs/source/gprmodelling.rst b/docs/source/gprmodelling.rst new file mode 100644 index 00000000..fecf2126 --- /dev/null +++ b/docs/source/gprmodelling.rst @@ -0,0 +1,114 @@ +.. _guidance: + +************************* +Guidance on GPR Modelling +************************* + +Basic concepts +============== + +This section discusses some basic concepts of GPR modelling. gprMax solves Maxwell's equations using the *Finite-Difference Time-Domain (FDTD)* method. Detailed information about the use of FDTD for GPR modelling can be found in [GIA1997]_ and general information on the FDTD method in [KUN1993]_ and [TAF2005]_. + +All electromagnetic phenomena, on a macroscopic scale, are described by the well-known Maxwell's equations. These are first order partial differential equations which express the relations between the fundamental electromagnetic field quantities and their dependence on their sources. + +.. math:: + + &\boldsymbol{\nabla}\boldsymbol{\times}\mathbf{E} =- \frac{\partial \mathbf{B}}{\partial t} \\ + &\boldsymbol{\nabla}\boldsymbol{\times}\mathbf{H} = \frac{\partial \mathbf{D}}{\partial t}+\mathbf{J_c}+\mathbf{J_s} \\ + &\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbf{B} = 0 \\ + &\boldsymbol{\nabla}\boldsymbol{\cdot}\mathbf{D} = q_v + +where :math:`t` is time (seconds) and :math:`q_v` is the volume electric charge density (coulombs/cubic metre). In Maxwell's equations, the field vectors are assumed to be single-valued, bounded, continuous functions of position and time. In order to simulate the GPR response from a particular target or set of targets the above equations have to be solved subject to the geometry of the problem and the initial conditions. + +The nature of the GPR forward problem classifies it as an *initial value -- open boundary* problem. This means that in order to obtain a solution you have to define an initial condition (i.e. excitation of the GPR transmitting antenna) and allow for the resulting fields to propagate through space reaching a zero value at infinity since, there is no specific boundary which limits the problem's geometry and where the electromagnetic fields can take a predetermined value. Although the first part is easy to accommodate (i.e. specification of the source), the second part can notbe easily tackled using a finite computational space. + +The FDTD approach to the numerical solution of Maxwell's equations is to discretize both the space and time continua. Thus the discretization spatial :math:`\Delta x`, :math:`\Delta y` and :math:`\Delta z` and +temporal :math:`\Delta t` steps play a very significant role -- since the smaller they are the closer the FDTD model is to a real representation of the problem. However, the values of the discretization steps always have to be finite, since computers have a limited amount of storage and finite processing speed. Hence, the FDTD model represents a discretized version of the real problem and is of limited size. The building block of this discretized FDTD grid is the Yee cell [YEE1966]_ named after Kane Yee who pioneered the FDTD method. This is illustrated for the 3D case in :numref:`yeecell`. + +.. _yeecell: + +.. figure:: images/yeecell.pdf + + 3D FDTD Yee cell + +By assigning appropriate constitutive parameters to the locations of the electromagnetic field components complex shaped targets can be included easily in the models. However, objects with curved boundaries are represented using a staircase approximation. + +The numerical solution is obtained directly in the time domain by using a discretized version of Maxwell's curl equations which are applied in each FDTD cell. Since these equations are discretized in both space and time the solution is obtained in an iterative fashion. In each iteration the electromagnetic fields advance (propagate) in the FDTD grid and each iteration corresponds to an elapsed simulated time of one :math:`\Delta t`. Hence by specifying the number of iterations you can instruct the FDTD solver to simulate the fields for a given time window. + +The price you has to pay of obtaining a solution directly in the time domain using the FDTD method is that the values of :math:`\Delta x`, :math:`\Delta y`, :math:`\Delta z` and :math:`\Delta t` can not be assigned independently. FDTD is a conditionally stable numerical process. The stability condition is known as the CFL condition after the initials of Courant, Freidrichs and Lewy and is given by, + +.. math:: \Delta t \leq \frac{1}{c\sqrt{\frac{1}{(\Delta x)^2}+\frac{1}{(\Delta y)^2}+\frac{1}{(\Delta z)^2}}}, + +where :math:`c` is the speed of light. Hence :math:`\Delta t` is bounded by the values of :math:`\Delta x`, :math:`\Delta y` and :math:`\Delta z`. The stability condition for the 2D case is easily obtained by letting :math:`\Delta z \longrightarrow \infty`. + +One of the most challenging issues in modelling *open boundary* problems, such as GPR, is the truncation of the computational domain at a finite distance from sources and targets where the values of the electromagnetic fields can not be calculated directly by the numerical method applied inside the model. Hence, an approximate condition known as *absorbing boundary condition (ABC)* is applied at a sufficient distance from the source to truncate and therefore limit the computational space. The role of this ABC is to absorb any waves impinging on it, hence simulating an unbounded space. The computational space (i.e the model) limited by the ABCs should contain all important features of the model such as sources and output points and targets. :numref:`abcs` illustrates this basic difference between the problem to be modelled and the actual FDTD modelled space. + +.. _abcs: + +.. figure:: images/abcs.pdf + + GPR forward problem showing computational domain bounded by Absorbing Boundary Conditions (ABCs) + +It is assumed that the half-space which contains the target(s) is of infinite extent. Therefore, the only reflected waves will be the ones originating from the target. In cases where the host medium is not of infinite extent (e.g. a finite concrete slab) the assumption of infinite extent can be made as far as the actual reflections from the slab termination are not of interest or its actual size is large enough that any reflected waves which will originate at its termination will not affect the solution for the required time window. In general, any objects that span the size of the computational domain (i.e. model) are assumed to extend to infinity. The only reflections which will originate from their termination at the truncation boundaries of the model are due to imperfections of the ABCs and in general are of a very small amplitude compared with the reflections from target(s) inside the model. + +All other *boundary conditions* which apply at interfaces between different media in the FDTD model are automatically enforced in gprMax. + +In order to make the most of gprMax in modelling GPR responses ultimately you should be familiar with the FDTD method on which these programs are based. There is a very large amount of information available in the relevant literature. Good starting points are [KUN1993]_ and [TAF2005]_ where as the specific application of FDTD to the GPR forward problem is described in [GIA1997]_. + + +Coordinate system and conventions +================================= + +A right-handed Cartesian coordinate system is used with the origin of space coordinates in the *lower left corner* at (0,0,0). :numref:`coord3d` illustrates the coordinate system of gprMax. Only one row of cells in the x direction is depicted. The space coordinates range from the left edge of the first cell to the right edge of the last one. Assuming that :math:`\Delta x = 1` metre, if you wanted to allocate a rectangle with its x dimension equal to 3 metres and its lower x coordinate at 1 then the x range would be [1..4]. The 3D cells allocated by gprMax would be [1..3]. In the 3D FDTD cell there are no field components located at the centre of the cell. Electric field components are tangential to, and magnetic field components normal to the interfaces between cells. The field components depicted in :numref:`coord3d` correspond to space coordinate 1. Source and output points defined in space coordinates are directly converted to cell coordinates and the corresponding field components. + +.. _coord3d: + +.. figure:: images/coord3d.pdf + + gprMax coordinate system and conventions. + +The actual positions of field components for a given set of space coordinates (x, y, z) are: + +.. math:: + + &E_x~(x+\frac{\Delta x}{2}, y, z) \\ + &E_y~(x, y+\frac{\Delta y}{2}, z) \\ + &E_z~(x, y, z+\frac{\Delta z}{2}) \\ + &H_x~(x, y+\frac{\Delta y}{2}, z+\frac{\Delta z}{2}) \\ + &H_y~(x+\frac{\Delta x}{2}, y, z+\frac{\Delta z}{2}) \\ + &H_z~(x+\frac{\Delta x}{2}, y+\frac{\Delta y}{2}, z) + +Hertzian dipole sources as well as other electric field excitations (i.e. voltage sources, transmission lines) are located at the corresponding electric field components. + + +Discretisation +============== + +There is no specific guideline for choosing the right discretization for a given problem. In general, it depends on the required accuracy, the frequency content of the source pulse and the size of the targets. Obviously, all targets present in a model must be adequately resolved. This means, for example, that a cylinder with radius equal to one or two spatial steps does not really look like a cylinder! + +An other important factor which influences the discretization is the errors associated with numerical induced dispersion. This means that contrary to the real world where electromagnetic waves propagate with the same velocity irrespectively of their direction and frequency (assuming no dispersive media and far-field conditions) in the discrete one this is not the case. This error (details can be found in [GIA1997]_ and [KUN1993]_) can be kept in a minimum if the following *rule-of-thumb* is satisfied: + +**The discretization step should be at least ten times smaller than the smallest wavelength of the propagating electromagnetic fields.** + +.. math:: \Delta l = \frac{\lambda}{10} + +Note that in general low-loss media wavelengths are much smaller compared to free space. + + +Absorbing boundary conditions +============================= + +The absorbing boundary conditions (ABCs) employed in gprMax will, in general, perform well (i.e. without introducing significant artificial reflections) if all sources and targets are kept at least 15 cells away from them. gprMax uses Perfectly Matched Layer (PML) ABCs based on a recursive integration approach to the complex frequency shifted (CFS) PML [GIA2012]_. A general formulation of this RIPML, which can be used to develop any order of PML, has been used to implement first and second order CFS stretching functions. One of the attractions of the RIPML is that it is easily applied as a correction to the field quantities after the complete FDTD grid has been updated using the standard FDTD update equations. + +The cells of the RIPML, which have a user adjustable thickness, very efficiently absorb most waves that propagate in them. Although, source and output points can be specified inside these cells **it is wrong to do so** from the point of view of correct modelling. The fields inside these cells are not of interest to GPR modelling. Placing sources inside these cells could have effects that have not been studied and will certainly provide erroneous results from a GPR modeller's point of view. The requirement to keep sources and targets at least 15 cells away for the PML has to be taken into account when deciding the size of the model domain. Additionally, free space (i.e. air) should be always included above a source for at least 15-20 cells in GPR models. Obviously, the more cells there are between observation points, sources, targets and the absorbing boundaries, the better the results will be. + +gprMax now offers the ability (for advanced users) to customise the parameters of the PML which allows its performance to be better optimised for specific applications. For further details see the :ref:`PML commands section `. + +This user guide, can not serve as an in depth tutorial and a review of the FDTD method. However, some useful hints and tips are given here in order to cover the most fundamental aspects of using an FDTD based program and avoid the most common errors. + +**References** + +.. [GIA1997] Giannopoulos, A. (1997). The investigation of Transmission-Line Matrix and Finite-Difference Time-Domain Methods for the Forward Problem of Ground Probing Radar, D.Phil thesis, Department of Electronics, University of York, UK +.. [KUN1993] Kunz, K. S., & Luebbers, R. J. (1993). The finite difference time domain method for electromagnetics. CRC press. +.. [TAF2005] Taflove, A., & Hagness, S. C. (2005). Computational electrodynamics. Artech house. +.. [YEE1966] Yee, K. S. (1966). Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media. IEEE Trans. Antennas Propag, 14(3), 302-307. diff --git a/docs/source/images/GSSI_1500_cylinder.png b/docs/source/images/GSSI_1500_cylinder.png new file mode 100644 index 00000000..3dc44421 Binary files /dev/null and b/docs/source/images/GSSI_1500_cylinder.png differ diff --git a/docs/source/images/GSSI_1500_cylinder_Bscan_results.pdf b/docs/source/images/GSSI_1500_cylinder_Bscan_results.pdf new file mode 100644 index 00000000..cb95440f Binary files /dev/null and b/docs/source/images/GSSI_1500_cylinder_Bscan_results.pdf differ diff --git a/docs/source/images/abcs.pdf b/docs/source/images/abcs.pdf new file mode 100644 index 00000000..ec1313b8 Binary files /dev/null and b/docs/source/images/abcs.pdf differ diff --git a/docs/source/images/antenna_MALA_1200.png b/docs/source/images/antenna_MALA_1200.png new file mode 100644 index 00000000..40e6b1ae Binary files /dev/null and b/docs/source/images/antenna_MALA_1200.png differ diff --git a/docs/source/images/antenna_MALA_1200_fs_results.pdf b/docs/source/images/antenna_MALA_1200_fs_results.pdf new file mode 100644 index 00000000..6cf64c02 Binary files /dev/null and b/docs/source/images/antenna_MALA_1200_fs_results.pdf differ diff --git a/docs/source/images/contsine.pdf b/docs/source/images/contsine.pdf new file mode 100644 index 00000000..c6c0b8d1 Binary files /dev/null and b/docs/source/images/contsine.pdf differ diff --git a/docs/source/images/coord3d.pdf b/docs/source/images/coord3d.pdf new file mode 100644 index 00000000..558573f1 Binary files /dev/null and b/docs/source/images/coord3d.pdf differ diff --git a/docs/source/images/cylinder_Ascan.pdf b/docs/source/images/cylinder_Ascan.pdf new file mode 100644 index 00000000..44cc6606 Binary files /dev/null and b/docs/source/images/cylinder_Ascan.pdf differ diff --git a/docs/source/images/cylinder_Bscan_results.pdf b/docs/source/images/cylinder_Bscan_results.pdf new file mode 100644 index 00000000..8913131e Binary files /dev/null and b/docs/source/images/cylinder_Bscan_results.pdf differ diff --git a/docs/source/images/cylinder_half_space_geo.pdf b/docs/source/images/cylinder_half_space_geo.pdf new file mode 100644 index 00000000..3bdd4fe5 Binary files /dev/null and b/docs/source/images/cylinder_half_space_geo.pdf differ diff --git a/docs/source/images/editables/abcs.pages b/docs/source/images/editables/abcs.pages new file mode 100644 index 00000000..ca9c9242 Binary files /dev/null and b/docs/source/images/editables/abcs.pages differ diff --git a/docs/source/images/editables/cylinder_half_space_geo.pages b/docs/source/images/editables/cylinder_half_space_geo.pages new file mode 100644 index 00000000..329f6a4f Binary files /dev/null and b/docs/source/images/editables/cylinder_half_space_geo.pages differ diff --git a/docs/source/images/gaussian.pdf b/docs/source/images/gaussian.pdf new file mode 100644 index 00000000..04cb9999 Binary files /dev/null and b/docs/source/images/gaussian.pdf differ diff --git a/docs/source/images/gaussiandot.pdf b/docs/source/images/gaussiandot.pdf new file mode 100644 index 00000000..78d33998 Binary files /dev/null and b/docs/source/images/gaussiandot.pdf differ diff --git a/docs/source/images/gaussiandotdot.pdf b/docs/source/images/gaussiandotdot.pdf new file mode 100644 index 00000000..53d7dd84 Binary files /dev/null and b/docs/source/images/gaussiandotdot.pdf differ diff --git a/docs/source/images/gaussiandotdotdot.pdf b/docs/source/images/gaussiandotdotdot.pdf new file mode 100644 index 00000000..6aae0197 Binary files /dev/null and b/docs/source/images/gaussiandotdotdot.pdf differ diff --git a/docs/source/images/gaussiandotdotnorm.pdf b/docs/source/images/gaussiandotdotnorm.pdf new file mode 100644 index 00000000..445c40e6 Binary files /dev/null and b/docs/source/images/gaussiandotdotnorm.pdf differ diff --git a/docs/source/images/gaussiandotnorm.pdf b/docs/source/images/gaussiandotnorm.pdf new file mode 100644 index 00000000..e677ea6e Binary files /dev/null and b/docs/source/images/gaussiandotnorm.pdf differ diff --git a/docs/source/images/gprMax_FB_logo.png b/docs/source/images/gprMax_FB_logo.png new file mode 100644 index 00000000..bf5f5d31 Binary files /dev/null and b/docs/source/images/gprMax_FB_logo.png differ diff --git a/docs/source/images/gprMax_logo.png b/docs/source/images/gprMax_logo.png new file mode 100644 index 00000000..8741166e Binary files /dev/null and b/docs/source/images/gprMax_logo.png differ diff --git a/docs/source/images/heterogeneous_soil.png b/docs/source/images/heterogeneous_soil.png new file mode 100644 index 00000000..b5639b46 Binary files /dev/null and b/docs/source/images/heterogeneous_soil.png differ diff --git a/docs/source/images/paraview_pipeline.png b/docs/source/images/paraview_pipeline.png new file mode 100644 index 00000000..9b833ea2 Binary files /dev/null and b/docs/source/images/paraview_pipeline.png differ diff --git a/docs/source/images/paraview_toolbar.png b/docs/source/images/paraview_toolbar.png new file mode 100644 index 00000000..4221cabd Binary files /dev/null and b/docs/source/images/paraview_toolbar.png differ diff --git a/docs/source/images/ricker.pdf b/docs/source/images/ricker.pdf new file mode 100644 index 00000000..5757e53b Binary files /dev/null and b/docs/source/images/ricker.pdf differ diff --git a/docs/source/images/sine.pdf b/docs/source/images/sine.pdf new file mode 100644 index 00000000..ba0dad5a Binary files /dev/null and b/docs/source/images/sine.pdf differ diff --git a/docs/source/images/yeecell.pdf b/docs/source/images/yeecell.pdf new file mode 100644 index 00000000..65a7225b Binary files /dev/null and b/docs/source/images/yeecell.pdf differ diff --git a/docs/source/includereadme.rst b/docs/source/includereadme.rst new file mode 100644 index 00000000..bb0877ff --- /dev/null +++ b/docs/source/includereadme.rst @@ -0,0 +1 @@ +.. include:: ../../gprMax/README.rst \ No newline at end of file diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100644 index 00000000..aa54b8fe --- /dev/null +++ b/docs/source/index.rst @@ -0,0 +1,38 @@ +***************** +gprMax User Guide +***************** + +.. toctree:: + :maxdepth: 2 + :caption: Introduction + + includereadme + features + gprmodelling + + +.. toctree:: + :maxdepth: 2 + :caption: Using gprMax + + input + geometry_snapshots + output + openmp_mpi + + +.. toctree:: + :maxdepth: 2 + :caption: Help and Support + + examples_2D + examples_3D + faqs + + +.. toctree:: + :maxdepth: 2 + :caption: Appendices + + app_waveforms + app_source \ No newline at end of file diff --git a/docs/source/input.rst b/docs/source/input.rst new file mode 100644 index 00000000..8e0c56a4 --- /dev/null +++ b/docs/source/input.rst @@ -0,0 +1,808 @@ +.. _commands: + +******************* +Input file commands +******************* + +An input file has to be supplied to gprMax which should contain all the necessary information to run a GPR model. The input file is an ASCII text file which can be prepared with any text editor or word-processing program. In the input file the hash character (``#``) is reserved and is used to denote the beginning of a command which will be passed to gprMax. The general syntax of commands is: + +.. code-block:: none + + #command_name: parameter1 parameter2 parameter3 ... + +A command and associated parameters should occupy a single line of the input file, and only one command per line is allowed. Hence, the first character of a line containing a command **must** be the hash character (``#``). If the line starts with **any other character** it is ignored by the program. Therefore, user comments or descriptions can be included in the input file. If a line starts with a hash character (``#``) the program will expect a valid command. If the name of the command is not correct the program will abandon execution and issue an error message. When a command requires more than one parameter then these should be separated using a white space character. + +The order of commands in the input file is not important with the exception of object construction commands. + +To describe the commands that can be used in the input file and their parameters the following conventions are used: + +* ``f`` means a real number which can be entered using either a ``[.]`` separating the integral from the decimal part, e.g. 1.5, or in scientific notation, e.g. 15e-1 or 0.15e1. +* ``i`` means an integer number. +* ``c`` means a single character, e.g. ``y``. +* ``str`` means a string of characters with **no** white spaces in between, e.g ``sand``. +* ``file`` means a filename. +* ``[ ]`` square brackets are used to indicate optional parameters. + +Unless otherwise specified, the SI system of units is used throughout gprMax: + +* All parameters associated with simulated space (i.e. size of model, spatial increments, etc...) should be specified in **metres**. +* All parameters associated with time (i.e. total simulation time, time instants, etc...) should be specified in **seconds**. +* All parameters denoting frequency should be specified in **Hertz**. +* All parameters associated with spatial coordinates in the model should be specified in **metres**. The origin of the coordinate system **(0,0)** is at the lower left corner of the model. + +It is important to note that gprMax converts spatial and temporal parameters given in **metres** and **seconds** respectively to integer values corresponding to **FDTD cell coordinates** and **iteration number**. Therefore, rounding to the nearest integer number of the user defined values is performed. + +The fundamental spatial and temporal discretization steps are denoted as :math:`\Delta x` , :math:`\Delta y`, :math:`\Delta z` and :math:`\Delta t` respectively. + +The commands have been grouped into six categories: + +* **Essential** - required to run any model, such as the domain size and spatial discretization +* **General** - provide further control over the model +* **Media** - used to introduce different materials into the model +* **Object construction** - used to build geometric shapes with different constitutive parameters +* **Excitation and output** - used to place source and output points in the model +* **PML** - provide advanced customisation and optimisation of the absorbing boundary conditions + +Essential commands +================== + +Most of the commands are optional but there are some essential commands which are necessary in order to construct any model. For example, none of the media and object commands are necessary to run a model. However, without specifying any objects in the model gprMax will simulate free space (air), which on its own, is not particularly useful for GPR modelling. If you have not specified a command which is essential in order to run a model, for example the size of the model, gprMax will terminate execution and issue an appropriate error message. + +The essential commands are: + +#domain: +-------- + +Allows you to specify the size of the model. The syntax of the command is: + +.. code-block:: none + + #domain: f1 f2 f3 + +where ``f1 f2 f3`` are the size of the model in the x, y, and z directions respectively. For example to specify a 500 x 500 x 1000mm model use: ``#domain: 0.5 0.5 1.0`` + +#dx_dy_dz: +---------- + +Allows you to specify the discretization of space in the x , y and z directions respectively (i.e. :math:`\Delta x` , :math:`\Delta y`, :math:`\Delta z`). The syntax of the command is: + +.. code-block:: none + + #dx_dy_dz: f1 f2 f3 + +where ``f1`` is the spatial step in the x direction (:math:`\Delta x`), ``f2`` is the spatial step in the y direction (:math:`\Delta y`) and ``f3`` is the spatial step in the z direction (:math:`\Delta z`). The spatial discretization controls the maximum permissible time step :math:`\Delta t` with which the solution advances in time in order to reach the required simulated time window. The relation between :math:`\Delta t` and :math:`\Delta x` , :math:`\Delta y`, :math:`\Delta z` is: + +.. math:: \Delta t \leq \frac{1}{c\sqrt{\frac{1}{(\Delta x)^2}+\frac{1}{(\Delta y)^2}+\frac{1}{(\Delta z)^2}}}, + +where :math:`c` is the speed of light. In gprMax the equality is used to determine :math:`\Delta t` from :math:`\Delta x` , :math:`\Delta y`, and :math:`\Delta z`. Small values of :math:`\Delta x` , :math:`\Delta y`, and :math:`\Delta z` result in small values for :math:`\Delta t` which means more iterations in order to reach a given simulated time. However, it is important to note that the smaller the values of :math:`\Delta x` , :math:`\Delta y`, :math:`\Delta z` and :math:`\Delta t` are the more accurate your model will be. See the :ref:`guidance` section for tips on choosing a spatial discretisation. + +#time_window: +------------- + +Allows you to specify the total required simulated time. The syntax of the command is: + +.. code-block:: none + + #time_window: f1 + +or + +.. code-block:: none + + #time_window: i1 + +In the first case the ``f1`` parameter determines the required simulated time in seconds. For example, if you want to simulate a GPR trace of 20 nanoseconds then ``#time_window: 20e-9`` can be used. gprMax will perform the necessary number of iterations in order to reach the required simulated time. Alternatively, if the command is specified with an ``i1`` gprMax will interpret this value as a total number of iterations. Hence the command ``#time_window: 100`` means that 100 iterations will be performed. The number of iterations and the total simulated time window are related by: + +.. math:: t_w = \Delta t × N_{it}, + +where :math:`t_w` is the time window in seconds, :math:`\Delta t` the time step, and :math:`N_{it}` the number of iterations. gprMax converts the specified time window in seconds to a number of iterations internally using the aforementioned equation. The result of the division is rounded to the nearest integer. + + +General commands +================ + +.. _python: + +#python: and #end_python: +------------------------- + +Allows you to write blocks of Python code between ``#python`` and ``#end_python`` in the input file. The code is executed when the input file is read by gprMax. + +For example, to use Python to automatically generate repetitive geometry: + +.. code-block:: none + + #python: + for x in range(0, 8) + print(’#cylinder: z 0.000 0.100 {} 0.050 0.005 pec’.format(0.020 + x * 0.020)) + #end_python: + +You can access the following built-in constants from your Python code: + +* ``c`` which is the speed of light in vacuum :math:`c=2.9979245 \times 10^8` m/s +* ``e0`` which is the permittivity of free space :math:`\epsilon_0=8.854187 \times 10^{-12}` F/m +* ``m0`` which is the permeability of free space :math:`\mu_0=1.256637 \times 10^{-6}` H/m +* ``z0`` which is the impedance of free space :math:`z_0=376.7303134` Ohms + +You can access the following built-in variables from your Python code: + +* ``current_model_run`` which is the current run number of the model that is been executed. +* ``number_model_runs`` which is the total number of runs specified when the model was initiatially executed. + +For example, if you running a model multiple times, e.g. to generate a scan, you can use the syntax: ``python -m gprMax my_input_file -n number_of_model_runs`` + +You can also access the built-in library of antenna models. Currently models of antennas similar to Geophysical Survey Systems, Inc. (GSSI) (http://www.geophysical.com) 1.5 GHz (Model 5100) antenna, and MALA Geoscience (http://www.malags.com/) 1.2 GHz antenna are included. They can be accessed from within a block of Python code using: ``from user_libs.antennas import antenna_like_GSSI_1500`` or ``from user_libs.antennas import antenna_like_MALA_1200``. + +For example, to use Python to include one of the antenna models from the built-in library at a location 0.125m, 0.094m, 0.100m (x,y,z) using a 1mm spatial resolution: + +.. code-block:: none + + #python: + from user_libs.antennas import antenna_like_GSSI_1500 + antenna_like_GSSI_1500(0.125, 0.094, 0.100, 0.001) + #end_python: + +#time_step_limit_type: +---------------------------- + +Allows you to choose whether to set the time step based on the 2D or 3D CFL condition. This command is useful when running a 2D model, i.e. when one dimension of the domain is one cell thick, because it allows the time step to be relaxed to the 2D CFL condition. The syntax of the command is: + +.. code-block:: none + + #time_step_limit_type: str1 + +where ``str1`` can be either 2D or 3D. + +#time_step_stability_factor: +---------------------------- + +Allows you to alter the value of the time step :math:`\Delta t` used by gprMax. gprMax uses the equality in the CFL condition, hence the maximum permissible time step. If a smaller time step is required then the syntax of the command is: + +.. code-block:: none + + #time_step_stability_factor: f1 + +where ``f1`` can take values :math:`0 < \textrm{f1} \leq 1`. Then the actual time step used will be :math:`\textrm{f1} \times \Delta t`, where :math:`\Delta t` is calculated using the equality from the CFL condition. + +#title: +------- + +Allows you to include a title for your model. This title is saved in the output file(s). The syntax of the command is: + +.. code-block:: none + + #title: str1 + +where ``str1`` can contain white space characters to separate individual words. The title has to be contained in a single line. + +#messages: +---------- + +Allows you to control the amount of information displayed on screen when gprMax is run. The syntax of the command is: + +.. code-block:: none + + #messages: c1 + +where ``c1`` can be either y (yes) or n (no) which turns on or off the messages on the screen. The default value is y. When messages are on, gprMax will display on the screen information the translation of space and time values to cell coordinates, iteration number, material parameters etc... This information can be useful for error checking. + +#num_threads: +----------------- + +Allows you to control how many OpenMP threads (usually the number of CPU cores available) are used when running the model. The most computationally intensive parts of gprMax, which are the FDTD solver loops, have been parallelised using OpenMP (http://openmp.org) which supports multi-platform shared memory multiprocessing. The syntax of the command is: + +.. code-block:: none + + #num_threads: i1 + +where ``i1`` is the number of OpenMP threads to use. If ``#num_threads`` is not specified gprMax will firstly look to see if the environment variable ``OMP_NUM_THREADS`` exists, and if not will detect and use all available CPU cores on the machine. + +.. _geometryview: + +#geometry_view: +--------------- + +Allows you output to file(s) information about the geometry of model. The file(s) use the open source Visualization ToolKit (VTK) (http://www.vtk.org) format which can be viewed in many free readers, such as Paraview (http://www.paraview.org). The command can be used to create several 3D views of the model which are useful for checking that it has been constructed as desired. The syntax of the command is: + +.. code-block:: none + + #geometry_view: f1 f2 f3 f4 f5 f6 f7 f8 f9 file1 c1 + +* ``f1 f2 f3`` are the lower left (x,y,z) coordinates of the volume of the geometry view in metres. +* ``f4 f5 f6`` are the upper right (x,y,z) coordinates of the volume of the geometry view in metres. +* ``f7 f8 f9`` are the spatial discretisation of the geometry view in metres. Typically these will be the same as the spatial discretisation of the model but they can be courser if desired. +* ``file1`` is the filename of the file where the geometry view will be stored. +* ``c1`` can be either n (normal) or f (fine) which specifies whether to output the geometry information on a per-cell basis (n) or a per-cell-edge basis (f). The fine mode should be reserved for viewing detailed parts of the geometry that occupy small volumes, as using this mode can generate geometry files with large file sizes. + +.. tip:: + + When you want to just check the geometry of your model, run gprMax using the optional command line argument ``--geometry-only``. This will build the model and produce any geometry view files, but will not run the simulation. + +#snapshot: +---------- + +Allows you to obtain information about the electromagnetic fields within a volume of the model at a given time instant. The file(s) use the open source Visualization ToolKit (VTK) (http://www.vtk.org) format which can be viewed in many free readers, such as Paraview (http://www.paraview.org). The syntax of this command is: + +.. code-block:: none + + #snapshot: f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 file1 + +or + +.. code-block:: none + + #snapshot: f1 f2 f3 f4 f5 f6 f7 f8 f9 i1 file1 + +* ``f1 f2 f3`` are the lower left (x,y,z) coordinates of the volume of the snapshot in metres. +* ``f4 f5 f6`` are the upper right (x,y,z) coordinates of the volume of the snapshot in metres. +* ``f7 f8 f9`` are the spatial discretisation of the snapshot in metres. +* ``f10`` or ``i1`` are the snapshot time in seconds (float number) or the iteration number (integer number) respectively which denote the point in time at which the snapshot will be taken. +* ``file1`` is the filename of the file where the snapshot will be stored. + +For example to save a snapshot of the electromagnetic fields in the model at a simulated time of 3 nanoseconds use: ``#snapshot: 0 0 0 1 1 1 0.1 0.1 0.1 3e-9 snap1`` + + +.. _materials: + +Media commands +============== + +Built-in materials +------------------ + +gprMax has two builtin materials which can be used by specifying the identifiers ``pec`` and ``free_space``. These simulate a perfect electric conductor and air, i.e. a non-magnetic material with :math:`\epsilon_r = 1`, :math:`\sigma = 0`, respectively. Additionally the identifiers ``grass`` and ``water`` are currently reserved for internal use and should not be used unless you intentionally want to change their properties. + +#material: +---------- + +Allows you to introduce a material into the model described by a set of constitutive parameters. The syntax of the command is: + +.. code-block:: none + + #material: f1 f2 f3 f4 str1 + +* ``f1`` is the relative permittivity, :math:`\epsilon_r` +* ``f2`` is the conductivity (Siemens/metre), :math:`\sigma` +* ``f3`` is the relative permeability, :math:`\mu_r` +* ``f4`` is the magnetic conductivity, :math:`\sigma_*` +* ``str1`` is an identifier for the material. + +For example ``#material: 3 0.01 1 0 my_sand`` creates a material called ``my_sand`` which has a relative permittivity (frequency independent) of :math:`\epsilon_r = 3`, a conductivity of :math:`\sigma = 0.01` S/m, and is non-magnetic, i.e. :math:`\mu_r = 1` and :math:`\sigma_* = 0` + + +#add_dispersion_debye: +---------------------- + +Allows you to add dispersive properties to an already defined ``#material`` based on a multiple pole Debye formulation (see :ref:`capabilities` section). The syntax of the command is: + +.. code-block:: none + + #add_dispersion_debye: i1 f1 f2 f3 f4 ... str1 + +* ``i1`` is the number of Debye poles. +* ``f1`` is the difference between the DC (static) relative permittivity and the relative permittivity at infinite frequency, i.e. :math:`\Delta \epsilon_r = \epsilon_{rs} - \epsilon_{r \infty}` , for the first Debye pole. +* ``f2`` is the relaxation time (seconds), :math:`\tau`, for the first Debye pole. +* ``f3`` is the difference between the DC (static) relative permittivity and the relative permittivity at infinite frequency, i.e. :math:`\Delta \epsilon_r = \epsilon_{rs} - \epsilon_{r \infty}` , for the second Debye pole. +* ``f4`` is the relaxation time (seconds), :math:`\tau`, for the second Debye pole. +* ``str1`` identifies the material to add the dispersive properties to. + +For example to create a model of water with a single Debye pole, :math:`\epsilon_{rs} = 80.1`, :math:`\epsilon_{r \inf} = 4.9` and :math:`\tau = 9.231\times 10^{-12}` seconds use: ``#material: 4.9 0.0 1.0 0.0 my_water`` and ``#add_dispersion_debye: 1 75.2 9.231e-12 my_water``. + +Important notes: + +* You can continue to add pairs of values for :math:`\Delta \epsilon_r` and :math:`\tau` for as many Debye poles as you have specified with ``i1``. +* The relative permittivity in the ``#material`` command should be given as the relative permittivity at infinite frequency, i.e. :math:`\epsilon_{r \infty}`. +* Values for the relaxation times should always be greater than the time step :math:`\Delta t` used in the model. gprMax checks and verifies this condition and if it does not then will exit raising an error. + + +#add_dispersion_lorenz: +----------------------- + +Allows you to add dispersive properties to an already defined ``#material`` based on a multiple pole Lorenz formulation (see :ref:`capabilities` section). The syntax of the command is: + +.. code-block:: none + + #add_dispersion_lorenz: i1 f1 f2 f3 f4 f5 f6 ... str1 + +* ``i1`` is the number of Lorenz poles. +* ``f1`` is the difference between the DC (static) relative permittivity and the relative permittivity at infinite frequency, i.e. :math:`\Delta \epsilon_r = \epsilon_{rs} - \epsilon_{r \infty}` , for the first Lorenz pole. +* ``f2`` is the relaxation time (seconds), :math:`\tau`, for the first Lorenz pole. +* ``f3`` is the damping factor for the first Lorenz pole. +* ``f4`` is the difference between the DC (static) relative permittivity and the relative permittivity at infinite frequency, i.e. :math:`\Delta \epsilon_r = \epsilon_{rs} - \epsilon_{r \infty}` , for the second Lorenz pole. +* ``f5`` is the relaxation time (seconds), :math:`\tau`, for the second Lorenz pole. +* ``f6`` is the damping factor for the second Lorenz pole. +* ``str1`` identifies the material to add the dispersive properties to. + +*For example...* + +Important notes: + +* You can continue to add triplets of values for :math:`\Delta \epsilon_r` and :math:`\tau` for as many Lorenz poles as you have specified with ``i1``. +* The relative permittivity in the ``#material`` command should be given as the relative permittivity at infinite frequency, i.e. :math:`\epsilon_{r \infty}`. +* Values for the relaxation times should always be greater than the time step :math:`\Delta t` used in the model. gprMax checks and verifies this condition and if it does not then will exit raising an error. + + +#add_dispersion_drude: +---------------------- + +Allows you to add dispersive properties to an already defined ``#material`` based on a multiple pole Drude formulation (see :ref:`capabilities` section). The syntax of the command is: + +.. code-block:: none + + #add_dispersion_drude: i1 f1 f2 f3 f4 f5 f6 ... str1 + +* ``i1`` is the number of Drude poles. +* ``f1`` is the difference between the DC (static) relative permittivity and the relative permittivity at infinite frequency, i.e. :math:`\Delta \epsilon_r = \epsilon_{rs} - \epsilon_{r \infty}` , for the first Drude pole. +* ``f2`` is the relaxation time (seconds), :math:`\tau`, for the first Drude pole. +* ``f3`` is the **x** for the first Drude pole. +* ``f4`` is the difference between the DC (static) relative permittivity and the relative permittivity at infinite frequency, i.e. :math:`\Delta \epsilon_r = \epsilon_{rs} - \epsilon_{r \infty}` , for the second Drude pole. +* ``f5`` is the relaxation time (seconds), :math:`\tau`, for the second Drude pole. +* ``f6`` is the **x** for the second Drude pole. +* ``str1`` identifies the material to add the dispersive properties to. + +*For example...* + +Important notes: + +* You can continue to add triplets of values for :math:`\Delta \epsilon_r` and :math:`\tau` for as many Drude poles as you have specified with ``i1``. +* The relative permittivity in the ``#material`` command should be given as the relative permittivity at infinite frequency, i.e. :math:`\epsilon_{r \infty}`. +* Values for the relaxation times should always be greater than the time step :math:`\Delta t` used in the model. gprMax checks and verifies this condition and if it does not then will exit raising an error. + + +#soil_peplinski: +---------------- + +Allows you to use a mixing model for soils proposed by Peplinski (http://dx.doi.org/10.1109/36.387598). The command is designed to be used in conjunction with the ``#fractal_box`` command for creating soils with realistic dielectric and geometric properties. The syntax of the command is: + +.. code-block:: none + + #soil_peplinski: f1 f2 f3 f4 f5 f6 str1 + +* ``f1`` is the sand fraction of the soil. +* ``f2`` is the clay fraction of the soil. +* ``f3`` is the bulk density of the soil in grams per centimetre cubed. +* ``f4`` is the density of the sand particles in the soil in grams per centimetre cubed. +* ``f5`` and ``f6`` define a range for the volumetric water fraction of the soil. +* ``str1`` is an identifier for the soil. + +For example for a soil with sand fraction 0.5, clay fraction 0.5, bulk density :math:`2~g/cm^3`, sand particle density of :math:`2.66~g/cm^3`, and a volumetric water fraction range of 0.001 - 0.25 use: ``#soil_peplinski: 0.5 0.5 2.0 2.66 0.001 0.25 my_soil``. + +Object construction commands +============================ + +Object construction commands are processed in the order they appear in the input file. Therefore space in the model allocated to a specific material using for example the ``#box`` command can be reallocated to another material using the same or any other object construction command. Space in the model can be regarded as a canvas in which objects are introduced and one can be overlaid on top of the other overwriting its properties in order to produce the desired geometry. The object construction commands can therefore be used to create complex shapes and configurations. + +Anisotropy +---------- + +It is possible to specify objects that have diagonal anisotropy which allows materials such as wood and fibre-reinforced composites, often imaged with GPR, to be more accurately modelled. + +.. math:: + + \bar{\bar{\epsilon}} = \left[ \begin{array}{ccc} + \epsilon_{xx} & 0 & 0 \\ + 0 & \epsilon_{yy} & 0 \\ + 0 & 0 & \epsilon_{zz} + \end{array} \right],\quad + \bar{\bar{\sigma}}= \left[ \begin{array}{ccc} + \sigma_{xx} & 0 & 0 \\ + 0 & \sigma_{yy} & 0 \\ + 0 & 0 & \sigma_{zz} + \end{array} \right] + +Standard isotropic objects specify one material identifier that defines the same properties in x, y, and z directions. However, every volumetric object building command can also be specified with three material identifiers, which allows properties for the x, y, and z directions to be separately defined. The ``#plate`` command, which defines a surface, can specify up to two material identifiers, and the ``#edge`` command, which defines a line, continues to take one material identifier. For example to create a box with different material properties in each of the x, y, and z directions use: + +.. code-block:: none + + #material: 41 10 1 0 matX + #material: 35 10 1 0 matY + #material: 33 1 1 0 matZ + #box: 0 0 0 0.1 0.1 0.1 matX matY matZ + +As another example, to create a cylinder of radius 10 mm that has the same properties in the x and y directions but different properties in the z direction use: + +.. code-block:: none + + #material: 41 10 1 0 matXY + #material: 33 1 1 0 matZ + #cylinder: 0.1 0.1 0.1 0.5 0.1 0.1 0.01 matXY matXY matZ + + +Dielectric smoothing +-------------------- + +At the boundaries between different materials in the model there is the question of which material properties to use. Should the last object to be defined at that location dictate the properties? Should an average set of properties of the materials of the objects that share that location be used? This latter option is often referred to as dielectric smoothing. To address this question gprMax includes an option to turn dielectric smoothing on or off for volumetric object building commands. The default behaviour (if no option is specified) is for dielectric smoothing to be on. The option can be specified with a single character ``y`` (on) or ``n`` (off) given after the material identifier in each object command. For example to specify a sphere of material ``sand`` with dielectric smoothing turned off use: ``#sphere: 0.5 0.5 0.5 0.1 sand n``. + +Important notes: + +* If an object is anistropic then dielectric smoothing is automatically turned off for that object. +* Non-volumetric object building commands, ``#edge`` and ``#plate`` cannot have dielectric smoothing. + +#edge: +------ + +Allows you to introduce a wire with specific properties into the model. A wire is an edge of a Yee cell and it can be useful to model resistors or thin wires. The syntax of the command is: + +.. code-block:: none + + #edge: f1 f2 f3 f4 f5 f6 str1 + +* ``f1 f2 f3`` are the starting (x,y,z) coordinates of the edge, and ``f4 f5 f6`` are the ending (x,y,z) coordinates of the edge. The coordinates should define a single line. +* ``str1`` is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials ``pec`` or ``free_space``. + +For example to specify a x-directed wire that is a perfect electric conductor, use: ``#edge: 0.5 0.5 0.5 0.7 0.5 0.5 pec``. Note that the y and z coordinates are identical. + +#plate: +------- + +Allows you to introduce a plate with specific properties into the model. A plate is a surface of a Yee cell and it can be useful to model objects thinner than a Yee cell. The syntax of the command is: + +.. code-block:: none + + #plate: f1 f2 f3 f4 f5 f6 str1 + +* ``f1 f2 f3`` are the lower left (x,y,z) coordinates of the plate, and ``f4 f5 f6`` are the upper right (x,y,z) coordinates of the plate. The coordinates should define a surface and not a 3D object like the ``#box`` command. +* ``str1`` is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials ``pec`` or ``free_space``. + +For example to specify a xy oriented plate that is a perfect electric conductor, use: ``#plate: 0.5 0.5 0.5 0.7 0.8 0.5 pec``. Note that the z coordinates are identical. + +#triangle: +---------- + +Allows you to introduce a triangular patch or a triangular prism with specific properties into the model. The patch is just a triangular surface made as a collection of staircased Yee cells, and the triangular prism extends the triangular patch in the direction perpendicular to the plane. The syntax of the command is: + +.. code-block:: none + + #triangle: f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 str1 + +* ``f1 f2 f3`` are the coordinates (x,y,z) of the first apex of the triangle, ``f4 f5 f6`` the coordinates (x,y,z) of the second apex, and ``f7 f8 f9`` the coordinates (x,y,z) of the third apex. +* ``f10`` is the thickness of the triangular prism. If the thickness is zero then a triangular patch is created. +* ``str1`` is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials ``pec`` or ``free_space``. + +For example, to specify a xy orientated triangular patch that is a perfect electric conductor, use: ``#triangle: 0.5 0.5 0.5 0.6 0.4 0.5 0.7 0.9 0.5 0.0 pec``. Note that the z coordinates are identical and the thickness is zero. + +#box: +----- + +Allows you to introduce an orthogonal parallelepiped with specific properties into the model. The syntax of the command is: + +.. code-block:: none + + #box: f1 f2 f3 f4 f5 f6 str1 + +* ``f1 f2 f3`` are the lower left (x,y,z) coordinates of the parallelepiped, and ``f4 f5 f6`` are the upper right (x,y,z) coordinates of the parallelepiped. +* ``str1`` is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials ``pec`` or ``free_space``. + +#sphere: +-------- + +Allows you to introduce a spherical object with specific parameters into the model. The syntax of the command is: + +.. code-block:: none + + #sphere: f1 f2 f3 f4 str1 + +* ``f1 f2 f3`` are the coordinates (x,y,z) of the centre of the sphere. +* ``f4`` is its radius. +* ``str1`` is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials ``pec`` or ``free_space``. + +For example, to specify a sphere with centre at (0.5, 0.5, 0.5), radius 100 mm, and with constitutive parameters of ``my_sand``, use: ``#sphere: 0.5 0.5 0.5 0.1 my_sand``. + +Important notes: + +* Sphere objects are permitted to extend outwith the model domain if desired, however, only parts of object inside the domain will be created. + +#cylinder: +---------- + +Allows you to introduce a circular cylinder into the model. The orientation of the cylinder axis can be arbitrary, i.e. it does not have align with one of the Cartesian axes of the model. The syntax of the command is: + +.. code-block:: none + + #cylinder: f1 f2 f3 f4 f5 f6 f7 str1 + +* ``f1 f2 f3`` are the coordinates (x,y,z) of the centre of one face of the cylinder, and ``f4 f5 f6`` are the coordinates (x,y,z) of the centre of the other face. +* ``f7`` is the radius of the cylinder. +* ``str1`` is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials ``pec`` or ``free_space``. + +For example, to specify a cylinder with its axis in the y direction, a length of 0.7 m, a radius of 100 mm, and that is a perfect electric conductor, use: ``#cylinder: 0.5 0.1 0.5 0.5 0.8 0.5 0.1 pec``. + +Important notes: + +* Cylinder objects are permitted to extend outwith the model domain if desired, however, only parts of object inside the domain will be created. + + +#cylindrical_sector: +-------------------- + +Allows you to introduce a cylindrical sector (shaped like a slice of pie) into the model. The syntax of the command is: + +.. code-block:: none + + #cylindrical_sector: c1 f1 f2 f3 f4 f5 f6 f7 str1 + +* ``c1`` is the direction of the axis of the cylinder from which the sector is defined and can be ``x``, ``y``, or ``z``. +* ``f1 f2`` are the coordinates of the centre of the cylindrical sector. +* ``f3 f4`` are the lower and higher coordinates of the axis of the cylinder from which the sector is defined (in effect they specify the thickness of the sector). +* ``f5`` is the radius of the cylindrical sector. +* ``f6`` is the starting angle (in degrees) for the cylindrical sector (with zero degrees defined on the positive first axis of the plane of the cylindrical sector). +* ``f7`` is the angle (in degrees) swept by the cylindrical sector (the finishing angle of the sector is always anti-clockwise from the starting angle). +* ``str1`` is a material identifier that must correspond to material that has already been defined in the input file, or is one of the builtin materials ``pec`` or ``free_space``. + +For example, to specify a cylindrical sector with its axis in the z direction, radius of 0.25 m, thickness of 2 mm, a starting angle of 330 :math:`^\circ`, a sector angle of 60 :math:`^\circ`, and that is a perfect electric conductor, use: ``#cylindrical_sector: z 0.34 0.24 0.500 0.502 0.25 330 60 pec``. + +Important notes: + +* Cylindrical sector objects are permitted to extend outwith the model domain if desired, however, only parts of object inside the domain will be created. + +.. _fractals: + +#fractal_box: +------------- + +Allows you to introduce an orthogonal parallelepiped with fractal distributed properties which are related to a mixing model or normal material into the model. The syntax of the command is: + +.. code-block:: none + + #fractal_box: f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 i1 str1 str2 [i2] + +* ``f1 f2 f3`` are the lower left (x,y,z) coordinates of the parallelepiped, and ``f4 f5 f6`` are the upper right (x,y,z) coordinates of the parallelepiped. +* ``f7`` is the fractal dimension which, for an orthogonal parallelepiped, should take values between zero and three. +* ``f8`` is used to weight the fractal in the x direction. +* ``f9`` is used to weight the fractal in the y direction. +* ``f10`` is used to weight the fractal in the z direction. +* ``i1`` is the number of materials to use for the fractal distribution (defined according to the associated mixing model). This should be set to one if using a normal material instead of a mixing model. +* ``str1`` is an identifier for the associated mixing model or material. +* ``str2`` is an identifier for the fractal box itself. +* ``i2`` is an optional parameter which controls the seeding of the random number generator used to create the fractals. By default (if you don't specify this parameter) the random number generator will be seeded by trying to read data from ``/dev/urandom`` (or the Windows analogue) if available or from the clock otherwise. + +For example, to create an orthogonal parallelepiped with fractal distributed properties using a Peplinski mixing model for soil, with 50 different materials over a range of water volumetric fractions from 0.001 - 0.25, you should first define the mixing model using: ``#soil_peplinski: 0.5 0.5 2.0 2.66 0.001 0.25 my_soil`` and then specify the fractal box using ``#fractal_box: 0 0 0 0.1 0.1 0.1 1.5 1 1 1 50 my_soil my_fractal_box``. + +#add_surface_roughness: +----------------------- + +Allows you to add rough surfaces to a ``#fractal_box`` in the model. A fractal distribution is used for the profile of the rough surface. The syntax of the command is: + +.. code-block:: none + + #add_surface_roughness: f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 str1 [i1] + +* ``f1 f2 f3`` are the lower left (x,y,z) coordinates of a surface on a ``#fractal_box``, and ``f4 f5 f6`` are the upper right (x,y,z) coordinates of a surface on a ``#fractal_box``. The coordinates must locate one of the six surfaces of a ``#fractal_box`` but do not have to extend over the entire surface. +* ``f7`` is the fractal dimension which, for an orthogonal parallelepiped, should take values between zero and three. +* ``f8`` is used to weight the fractal in the first direction of the surface. +* ``f9`` is used to weight the fractal in the second direction of the surface. +* ``f10 f11`` define lower and upper limits for a range over which the roughness can vary. These limits should be specified relative to the dimensions of the ``#fractal_box`` that the rough surface is being applied. +* ``str1`` is an identifier for the ``#fractal_box`` that the rough surface should be applied to. +* ``i1`` is an optional parameter which controls the seeding of the random number generator used to create the fractals. By default (if you don't specify this parameter) the random number generator will be seeded by trying to read data from ``/dev/urandom`` (or the Windows analogue) if available or from the clock otherwise. + +Up to six ``#add_rough_surface commands`` can be given for any ``#fractal_box`` corresponding to the six surfaces. + +For example, if a ``#fractal_box`` has been specified using: ``#fractal_box: 0 0 0 0.1 0.1 0.1 1.5 1 1 1 50 my_soil my_fractal_box`` then to apply a rough surface that varys between 85 mm and 110 mm (i.e. valleys that are up to 15 mm deep and peaks that are up to 10 mm tall) to the surface that is in the positive z direction, use ``#add_surface_roughness: 0 0 0.1 0.1 0.1 0.1 1.5 1 1 0.085 0.110 my_fractal_box``. + +#add_surface_water: +------------------- + +Allows you to add surface water to a ``#fractal_box`` in the model that has had a rough surface applied. The syntax of the command is: + +.. code-block:: none + + #add_surface_water: f1 f2 f3 f4 f5 f6 f7 str1 + +* ``f1 f2 f3`` are the lower left (x,y,z) coordinates of a surface on a ``#fractal_box``, and ``f4 f5 f6`` are the upper right (x,y,z) coordinates of a surface on a ``#fractal_box``. The coordinates must locate one of the six surfaces of a ``#fractal_box`` but do not have to extend over the entire surface. +* ``f7`` defines the depth of the water, which should be specified relative to the dimensions of the ``#fractal_box`` that the surface water is being applied. +* ``str1`` is an identifier for the ``#fractal_box`` that the surface water should be applied to. + +For example, to add surface water that is 5 mm deep to an existing ``#fractal_box`` that has been specified using ``#fractal_box: 0 0 0 0.1 0.1 0.1 1.5 1 1 1 50 my_soil my_fractal_box`` and has had a rough surface applied using ``#add_surface_roughness: 0 0 0.1 0.1 0.1 0.1 1.5 1 1 0.085 0.110 my_fractal_box``, use ``#add_surface_water: 0 0 0.1 0.1 0.1 0.1 0.105 my_fractal_box``. + +Important notes: + +* The water is modelled using a single-pole Debye formulation with properties :math:`\epsilon_{rs} = 80.1`, :math:`\epsilon_{\infty} = 4.9`, and a relaxation time of :math:`\tau = 9.231 \times 10^{-12}` seconds (http://dx.doi.org/10.1109/TGRS.2006.873208). If you prefer, gprMax will use your own definition for water as long as it is named ``water``. + +#add_grass: +----------- + +Allows you to add grass with roots to a ``#fractal_box`` in the model. The blades of grass are randomly distributed over the specified surface area and a fractal distribution is used to vary the height of the blades of grass and depth of the grass roots. The syntax of the command is: + +.. code-block:: none + + #add_grass: f1 f2 f3 f4 f5 f6 f7 f8 f9 i1 str1 [i2] + +* ``f1 f2 f3`` are the lower left (x,y,z) coordinates of a surface on a ``#fractal_box``, and ``f4 f5 f6`` are the upper right (x,y,z) coordinates of a surface on a ``#fractal_box``. The coordinates must locate one of three surfaces (in the positive axis direction) of a ``#fractal_box`` but do not have to extend over the entire surface. +* ``f7`` is the fractal dimension which, for an orthogonal parallelepiped, should take values between zero and three. +* ``f8 f9`` define lower and upper limits for a range over which the height of the blades of grass can vary. These limits should be specified relative to the dimensions of the ``#fractal_box`` that the grass is being applied. +* ``i1`` is the number of blades of grass that should be applied to the surface area. +* ``str1`` is an identifier for the ``#fractal_box`` that the grass should be applied to. +* ``i2`` is an optional parameter which controls the seeding of the random number generator used to create the fractals. By default (if you don't specify this parameter) the random number generator will be seeded by trying to read data from ``/dev/urandom`` (or the Windows analogue) if available or from the clock otherwise. + +For example, to apply 100 blades of grass that vary in height between 100 and 150 mm to the entire surface in the positive z direction of a ``#fractal_box`` that had been specified using ``#fractal_box: 0 0 0 0.1 0.1 0.1 1.5 1 1 50 my_soil my_fractal_box``, use: ``#add_grass: 0 0 0.1 0.1 0.1 0.1 1.5 0.2 0.25 100 my_fractal_box``. + +Important notes: + +* The grass is modelled using a single-pole Debye formulation with properties :math:`\epsilon_{rs} = 18.5087`, :math:`\epsilon_{\infty} = 12.7174`, and a relaxation time of :math:`\tau = 1.0793 \times 10^{-11}` seconds (http://dx.doi.org/10.1007/BF00902994). If you prefer, gprMax will use your own definition for grass if you use a material named ``grass``. The geometry of the blades of grass are defined by the parametric equations: :math:`x = x_c +s_x {\left( \frac{t}{b_x} \right)}^2`, :math:`y = y_c +s_y {\left( \frac{t}{b_y} \right)}^2`, and :math:`z=t`, where :math:`s_x` and :math:`s_y` can be -1 or 1 which are randomly chosen, and where the constants :math:`b_x` and :math:`b_y` are random numbers based on a Gaussian distribution. + +Excitation commands +=================== + +#waveform: +---------- + +Allows you to specify waveforms to use with sources in the model. The syntax of the command is: + +.. code-block:: none + + #waveform: str1 f1 f2 str2 + +* ``str1`` is the type of waveform which can be: + + * ``gaussian`` which is a Gaussian waveform. + * ``gaussiandot`` which is the first derivative of a Gaussian waveform. + * ``gaussiandotnorm`` which is the normalised first derivative of a Gaussian waveform. + * ``gaussiandotdot`` which is the second derivative of a Gaussian waveform. + * ``gaussiandotdotnorm`` which is the normalised second derivative of a Gaussian waveform. + * ``gaussiandotdotdot`` which is the third derivative of a Gaussian waveform. + * ``ricker`` which is a Ricker (or Mexican hat) waveform, i.e. the negative, normalised second derivative of a Gaussian waveform. + * ``sine`` which is a single cycle of a sine waveform. + * ``contsine`` which is a continuous sine waveform. In order to avoid introducing noise into the calculation the amplitude of the waveform is modulated for the first cycle of the sine wave (ramp excitation). +* ``f1`` is the amplitude of the waveform. +* ``f2`` is the frequency of the waveform in Hertz. +* ``str2`` is an identifier for the waveform used to assign it to a source. + +For example, to specify a Gaussian waveform with an amplitude of one and a centre frequency of 1.2 GHz, use: ``#waveform: gaussian 1 1.2e9 my_gauss_pulse``. + +Important notes: + +* The functions used to create the waveforms can be found in the :ref:`waveforms` appendix. + +#excitation_file: +----------------- + +Allows you to specify an ASCII file that contains columns of amplitude values that specify custom waveform shapes that can be used with sources in the model. The first row of each column must begin with a identifier string that will be used as the name of each waveform. There needs to be at least as many amplitude values as the number of iterations that are going to be performed. If there are less than this number then at the end of the sequence of amplitude values zero values will be added to pad the sequence up to the number of iterations. If extra amplitude values are specified than needed then they are ignored. The syntax of the command is: + +.. code-block:: none + + #excitation_file: file1 + +``str1`` is the name of the ASCII file containing the specified waveform. It should be located in the same directory as the input file. + +For example, to specify the file ``my_waves.txt``, which contains two custom waveform shapes, use: ``#excitation_file: my_waves.txt``. The contents of the file ``my_waves.txt`` would take the form: + +.. code-block:: none + + my_pulse1 my_pulse2 + 0 0 + 1.2e-6 0 + 1.3e-6 1e-1 + 5e-6 1.5e-1 + ... ... + ... ... + ... ... + +#hertzian_dipole: +----------------- + +Allows you to specify a current density term at an electric field location (the simplest excitation). This will simulate an infinitesimal electric dipole (it does have a length of :math:`\Delta l`). This is often referred to as an additive or soft source. The syntax of the command is: + +.. code-block:: none + + #hertzian_dipole: c1 f1 f2 f3 [f4 f5] str1 + +* ``c1`` is the polarisation of the source and can be ``x``, ``y``, or ``z``. +* ``f1 f2 f3`` are the coordinates (x,y,z) of the source in the model. +* ``f4 f5`` are optional parameters. ``f4`` is a time delay in starting the source. ``f5`` is a time to remove the source. If the time window is longer than the source removal time then the source will stop after the source removal time. If the source removal time is longer than the time window then the source will be active for the entire time window. If ``f4 f5`` are omitted the source will start at the beginning of time window and stop at the end of the time window. +* ``str1`` is the identifier of the waveform that should be used with the source. + +For example, to use a x polarised Hertzian dipole with unit amplitude and a 600 MHz centre frequency Ricker waveform, use: ``#waveform: ricker 1 600e6 my_ricker_pulse`` and ``#hertzian_dipole: x 0.05 0.05 0.05 my_ricker_pulse``. + +#magnetic_dipole: +----------------- + +This will simulate an infinitesimal magnetic dipole. This is often referred to as an additive or soft source. The syntax of the command is: + +.. code-block:: none + + #magnetic_dipole: c1 f1 f2 f3 [f4 f5] str1 + +* ``c1`` is the polarisation of the source and can be ``x``, ``y``, or ``z``. +* ``f1 f2 f3`` are the coordinates (x,y,z) of the source in the model. +* ``f4 f5`` are optional parameters. ``f4`` is a time delay in starting the source. ``f5`` is a time to remove the source. If the time window is longer than the source removal time then the source will stop after the source removal time. If the source removal time is longer than the time window then the source will be active for the entire time window. If ``f4 f5`` are omitted the source will start at the beginning of time window and stop at the end of the time window. +* ``str1`` is the identifier of the waveform that should be used with the source. + +#voltage_source: +---------------- + +Allows you to introduce a voltage source at the position of an electric field component either as a hard source (i.e. replacing the value of electic field component) or, as having an internal lumped resistance. It is useful for exciting GPR antennas when the physical properties of the antenna are included in the model. The syntax of the command is: + +.. code-block:: none + + #voltage_source: c1 f1 f2 f3 f4 [f5 f6] str1 + +* ``c1`` is the polarisation of the source and can be ``x``, ``y``, or ``z``. +* ``f1 f2 f3`` are the coordinates (x,y,z) of the source in the model. +* ``f4`` is the internal resistance of the voltage source in Ohms. If ``f4`` is set to zero then the voltage source is a hard source. That means it prescribes the value of the electric field component. If the waveform becomes zero then the source is perfectly reflecting. +* ``f5 f6`` are optional parameters. ``f5`` is a time delay in starting the source. ``f6`` is a time to remove the source. If the time window is longer than the source removal time then the source will stop after the source removal time. If the source removal time is longer than the time window then the source will be active for the entire time window. If ``f5 f6`` are omitted the source will start at the beginning of time window and stop at the end of the time window. +* ``str1`` is the identifier of the waveform that should be used with the source. + +For example, to specify a y directed voltage source with an internal resistance of 50 Ohms, an amplitude of five, and a 1.2 GHz centre frequency Gaussian waveform use: ``#waveform: gaussian 5 1.2e9 my_gauss_pulse`` and ``#voltage_source: y 0.05 0.05 0.05 50 my_gauss_pulse``. + +#rx: +---- + +Allows you to introduce output points into the model. These are locations where the values of the electric and magnetic field components over the number of iterations of the model will be saved to file. The syntax of the command is: + +.. code-block:: none + + #rx: f1 f2 f3 + +``f1 f2 f3`` are the coordinates (x,y,z) of the receiver in the model. + +#rx_box: +-------- + +Provides a simple method of defining multiple output points in the model. The syntax of the command is: + +.. code-block:: none + + #rx_box: f1 f2 f3 f4 f5 f6 f7 f8 f9 + +* ``f1 f2 f3`` are the lower left (x,y,z) coordinates of the output volume, and ``f4 f5 f6`` are the upper right (x,y,z) coordinates of the output volume. +* ``f7 f8 f9`` are the increments (x,y,z) which define the number of output points in each direction. The minimum value of ``f7`` is :math:`\Delta x`, the minimum value of ``f8`` is :math:`\Delta y`, and the minimum value of ``f9`` is :math:`\Delta z`. + +#src_steps: and #rx_steps: +-------------------------- + +Provide a simple method to allow you to move the location of all sources (``#src_steps``) or all receivers (``#rx_steps``) between runs of a model. The syntax of the commands is: + +.. code-block:: none + + #src_steps: f1 f2 f3 + #rx_steps: f1 f2 f3 + +``f1 f2 f3`` are increments (x,y,z) to move all sources (``#hertzian_dipole``, ``#magnetic_dipole``, or ``#voltage_source``) or all receivers (created using either ``#rx`` or ``#rx_box`` commands). + +.. _pml: + +PML commands +============ + +The default behaviour is for gprMax to use a first order CFS PML that has a thickness of 10 cells on each of the six sides of the model domain. This can be altered by using the following commands. + +#pml_cells: +------------ + +Allows you to control the number of cells of PML that are used on the six sides of the model domain. The PML is defined within the model domain, i.e. it is not added to the domain size. The syntax of the command is: + +.. code-block:: none + + #pml_cells: i1 [i2 i3 i4 i5 i6] + +* ``i1`` is the number of cells of PML to use on all sides of the model domain, or ``i1`` is the number of cells of PML to use on the side of the model domain in the negative x-axis direction. +* ``i2`` is the number of cells of PML to use on the side of the model domain in the negative y-axis direction. +* ``i3`` is the number of cells of PML to use on the side of the model domain in the negative z-axis direction. +* ``i4`` is the number of cells of PML to use on the side of the model domain in the positive x-axis direction. +* ``i5`` is the number of cells of PML to use on the side of the model domain in the positive y-axis direction. +* ``i6`` is the number of cells of PML to use on the side of the model domain in the positive z-axis direction. +* ``i1 i2 i3 i4 i5 i6`` may be set to zero to turn off the PML on a specific side of the model domain. + +To create a 2D model (a one cell slice of 3D) switch off the PML in the one cell (infinite) direction, e.g. to create a 2D model in the x-y plane, with a dimension of one cell in the z direction, and 10 cells of PML elsewhere, use the command: + +.. code-block:: none + + #pml_cells: 10 10 0 10 10 0 + +#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 f1 f2 str2 f3 f4 str3 f5 f6 + +* ``str1`` is the type of scaling to use for the CFS :math:`\alpha` parameter. It can be ``constant``, ``linear``, ``inverselinear``, ``quadratic``, ``cubic``, and ``quartic``. +* ``f1 f2`` are the minimum and maximum values for the CFS :math:`\alpha` parameter. +* ``str2`` is the type of scaling to use for the CFS :math:`\kappa` parameter. It can be ``constant``, ``linear``, ``inverselinear``, ``quadratic``, ``cubic``, and ``quartic``. +* ``f3 f4`` are the minimum and maximum values for the CFS :math:`\kappa` parameter. +* ``str3`` is the type of scaling to use for the CFS :math:`\sigma` parameter. It can be ``constant``, ``linear``, ``inverselinear``, ``quadratic``, ``cubic``, and ``quartic``. +* ``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 first order PML are: ``#pml_cfs: constant 0 0 constant 1 1 quartic 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. diff --git a/docs/source/models/GSSI_1500_cylinder.vti b/docs/source/models/GSSI_1500_cylinder.vti new file mode 100644 index 00000000..3d85e40d Binary files /dev/null and b/docs/source/models/GSSI_1500_cylinder.vti differ diff --git a/docs/source/models/GSSI_1500_cylinder_Bscan.in b/docs/source/models/GSSI_1500_cylinder_Bscan.in new file mode 100755 index 00000000..70b7862e --- /dev/null +++ b/docs/source/models/GSSI_1500_cylinder_Bscan.in @@ -0,0 +1,16 @@ +#title: B-scan from a metal cylinder buried in a dielectric half-space with a GSSI 1.5GHz 'like' antenna +#domain: 0.480 0.148 0.235 +#dx_dy_dz: 0.001 0.001 0.001 +#time_window: 6e-9 + +#material: 6 0 1 0 half_space + +#box: 0 0 0 0.480 0.148 0.170 half_space +#cylinder: 0.240 0 0.080 0.240 0.148 0.080 0.010 pec + +#python: +from user_libs.antennas import antenna_like_GSSI_1500 +antenna_like_GSSI_1500(0.105 + current_model_run * 0.005, 0.074, 0.170, 0.001) +#end_python: + +geometry_view: 0 0 0 0.480 0.148 0.235 0.001 0.001 0.001 GSSI_1500_cylinder n \ No newline at end of file diff --git a/docs/source/models/GSSI_1500_cylinder_Bscan_all.out b/docs/source/models/GSSI_1500_cylinder_Bscan_all.out new file mode 100644 index 00000000..de60bf61 Binary files /dev/null and b/docs/source/models/GSSI_1500_cylinder_Bscan_all.out differ diff --git a/docs/source/models/antenna_MALA_1200_fs.in b/docs/source/models/antenna_MALA_1200_fs.in new file mode 100755 index 00000000..a0187229 --- /dev/null +++ b/docs/source/models/antenna_MALA_1200_fs.in @@ -0,0 +1,9 @@ +#title: MALA 1.2GHz 'like' antenna in free-space +#domain: 0.264 0.189 0.220 +#dx_dy_dz: 0.001 0.001 0.001 +#time_window: 6e-9 + +#python: +from user_libs.antennas import antenna_like_MALA_1200 +antenna_like_MALA_1200(0.132, 0.095, 0.100, 0.001) +#end_python: \ No newline at end of file diff --git a/docs/source/models/antenna_MALA_1200_fs.out b/docs/source/models/antenna_MALA_1200_fs.out new file mode 100644 index 00000000..9badb9bb Binary files /dev/null and b/docs/source/models/antenna_MALA_1200_fs.out differ diff --git a/docs/source/models/cylinder_Ascan_2D.in b/docs/source/models/cylinder_Ascan_2D.in new file mode 100644 index 00000000..b2fa6c76 --- /dev/null +++ b/docs/source/models/cylinder_Ascan_2D.in @@ -0,0 +1,17 @@ +#title: A-scan from a metal cylinder buried in a dielectric half-space +#domain: 0.240 0.190 0.001 +#dx_dy_dz: 0.001 0.001 0.001 +#time_window: 4e-9 +#time_step_limit_type: 2D +#pml_cells: 10 10 0 10 10 0 + +#material: 6 0 1 0 half_space + +#waveform: ricker 1 1.5e9 my_ricker +#hertzian_dipole: z 0.100 0.170 0 my_ricker +#rx: 0.140 0.170 0 + +#box: 0 0 0 0.240 0.170 0.001 half_space +#cylinder: 0.120 0.080 0 0.120 0.080 0.001 0.010 pec + +#geometry_view: 0 0 0 0.240 0.190 0.001 0.001 0.001 0.001 cylinder_half_space n \ No newline at end of file diff --git a/docs/source/models/cylinder_Ascan_2D.out b/docs/source/models/cylinder_Ascan_2D.out new file mode 100644 index 00000000..0fcb075f Binary files /dev/null and b/docs/source/models/cylinder_Ascan_2D.out differ diff --git a/docs/source/models/cylinder_Bscan_2D.in b/docs/source/models/cylinder_Bscan_2D.in new file mode 100644 index 00000000..9609128d --- /dev/null +++ b/docs/source/models/cylinder_Bscan_2D.in @@ -0,0 +1,17 @@ +#title: B-scan from a metal cylinder buried in a dielectric half-space +#domain: 0.240 0.190 0.001 +#dx_dy_dz: 0.001 0.001 0.001 +#time_window: 4e-9 +#time_step_limit_type: 2D +#pml_cells: 10 10 0 10 10 0 + +#material: 6 0 1 0 half_space + +#waveform: ricker 1 1.5e9 my_ricker +#hertzian_dipole: z 0.020 0.170 0 my_ricker +#rx: 0.060 0.170 0 +#src_steps: 0.002 0 0 +#rx_steps: 0.002 0 0 + +#box: 0 0 0 0.240 0.170 0.001 half_space +#cylinder: 0.120 0.080 0 0.120 0.080 0.001 0.010 pec \ No newline at end of file diff --git a/docs/source/models/cylinder_Bscan_2D_all.out b/docs/source/models/cylinder_Bscan_2D_all.out new file mode 100644 index 00000000..26cba768 Binary files /dev/null and b/docs/source/models/cylinder_Bscan_2D_all.out differ diff --git a/docs/source/models/cylinder_half_space.vti b/docs/source/models/cylinder_half_space.vti new file mode 100644 index 00000000..08f7cbb1 Binary files /dev/null and b/docs/source/models/cylinder_half_space.vti differ diff --git a/docs/source/models/heterogeneous_soil.in b/docs/source/models/heterogeneous_soil.in new file mode 100755 index 00000000..4f0a10d5 --- /dev/null +++ b/docs/source/models/heterogeneous_soil.in @@ -0,0 +1,14 @@ +#title: Heterogeneous soil using a stochastic distribution of dielectric properties given by a mixing model from Peplinski +#domain: 0.15 0.15 0.1 +#dx_dy_dz: 0.001 0.001 0.001 +#time_window: 6e-9 + +#waveform: ricker 1 1.5e9 my_ricker +#hertzian_dipole: y 0.045 0.075 0.085 my_ricker +#rx: 0.105 0.075 0.085 + +#soil_peplinski: 0.5 0.5 2.0 2.66 0.001 0.25 my_soil +#fractal_box: 0 0 0 0.15 0.15 0.070 1.5 1 1 1 50 my_soil my_soil_box +#add_surface_roughness: 0 0 0.070 0.15 0.15 0.070 1.5 1 1 0.065 0.080 my_soil_box + +#geometry_view: 0 0 0 0.15 0.15 0.1 0.001 0.001 0.001 heterogeneous_soil n \ No newline at end of file diff --git a/docs/source/openmp_mpi.rst b/docs/source/openmp_mpi.rst new file mode 100644 index 00000000..0cdeec4c --- /dev/null +++ b/docs/source/openmp_mpi.rst @@ -0,0 +1,74 @@ +.. _openmp_mpi: + +************************ +Parallelism - OpenMP/MPI +************************ + +OpenMP +====== + +The most computationally intensive parts of gprMax, which are the FDTD solver loops, have been parallelised using OpenMP (http://openmp.org) which supports multi-platform shared memory multiprocessing. + +By default gprMax will try to lookup and use the maximum number of OpenMP threads (usually the number of CPU cores) available on your machine. You can override this behaviour in two ways: firstly, gprMax will check to see if the ``#num_threads`` command is present in your input file; if not, gprMax will check to see if the environment variable ``OMP_NUM_THREADS`` is set. This can be useful if you are running gprMax on a cluster or in a HPC environment where you might not want to use all of the available CPU cores. + +MPI +=== + +The Message Passing Interface (MPI) has been utilised to implement a simple task farm that can be used to distribute a series of models as independent tasks. This can be useful in many GPR simulations where a B-scan (composed of multiple A-scans) is required. Each A-scan can be task-farmed as a independent model. Within each independent model OpenMP threading will continue to be used. Overall this creates what is know as a mixed mode OpenMP/MPI job. + +By default the MPI task farm functionality is turned off. It can be switched on using the ``-mpi`` command line flag. MPI requires an installation of OpenMPI (http://www.open-mpi.org) and the mpi4py Python package. + +.. note:: + + It seems due to a lack of interest there are no binary versions of OpenMPI available for Microsoft Windows (https://www.open-mpi.org/software/ompi/v1.6/ms-windows.php) + +Running gprMax using the MPI task farm functionality is heavily dependent on the configuration of your machine/cluster. The following example is intended as general guidance to help you get started. + +Grid Engine example +------------------- + +Clusters usually requires jobs to be submitted to a queue using a job script. Typically within that script the ``mpirun`` program is used to execute MPI jobs. Here is an example of a job script for a cluster that uses Oracle (Sun) Grid Engine. The behaviour of most of the variables is explained in the comments in the script. + +.. code-block:: none + + #!/bin/bash + ##################################################################################### + ### Specify bash shell: + #$ -S /bin/bash + + ### Change to current working directory: + #$ -cwd + + ### Specify runtime (hh:mm:ss): + #$ -l h_rt=01:00:00 + + ### Email options: + #$ -m ea -M joe.bloggs@email.com + + ### Parallel environment ($NSLOTS): + #$ -pe openmpi_fillup_mark2 80 + + ### Job script name: + #$ -N test_mpi.sh + ##################################################################################### + + ### Initialise environment module + . /etc/profile.d/modules.sh + + ### Load Anaconda environment with Python 3 and packages + module load anaconda + source activate python3 + + ### Load OpenMPI + module load openmpi-gcc + + ### Set number of OpenMP threads + export OMP_NUM_THREADS=8 + + ### Run gprMax with input file + cd $HOME/gprMax + mpirun -np $NSLOTS python -m gprMax test.in -n 10 -mpi + +The ``NSLOTS`` variable is usually the number of MPI tasks multiplied by the number of OpenMP threads per task. In this example the number of MPI tasks is 10 and number of OpenMP threads per task is 8, so 8O slots are required. + + diff --git a/docs/source/output.rst b/docs/source/output.rst new file mode 100644 index 00000000..d1d0fe25 --- /dev/null +++ b/docs/source/output.rst @@ -0,0 +1,81 @@ +.. _output: + +*********** +Output file +*********** + +gprMax produces an output file that has the same name as the input file but with ``.out`` appended. The output file uses the widely-supported HDF5 (https://www.hdfgroup.org/HDF5/) format which was designed to store and organize large amounts of numerical data. + + +File structure +============== + +The output file has the following HDF5 attributes at the root (``/``): + +* ``Title`` is the title of the model +* ``Iterations`` is the number of iterations for the time window of the model +* ``dx, dy, dz`` is a tuple containing the spatial discretisation, i.e. :math:`\Delta x`, :math:`\Delta y`, :math:`\Delta z` +* ``dt`` is the time step of the model, i.e. :math:`\Delta t` +* ``srcsteps`` is the spatial increment used to move all sources between model runs. +* ``rxsteps`` is the spatial increment used to move all receivers between model runs. +* ``ntx`` is the total number of sources in the model. +* ``nrx`` is the total number of receievers in the model. + +The output file contains HDF5 groups for sources (``txs``) and receivers (``rxs``). Within each group are further groups that correspond to individual sources, e.g. ``tx1``, ``tx2`` etc..., and receivers, e.g. ``rx1``, ``rx2`` etc... + +.. code-block:: none + + / + rxs/ + rx1/ + Position + Ex + Ey + Ez + Hx + Hy + Hz + rx2/ + ... + txs/ + tx1/ + Position + tx2/ + ... + +Within each individual ``rx`` group are the following datasets: + +* ``Position`` is the x, y, z position (in metres) of the receiver in the model. +* ``Ex`` is an array containing the time history (for the model time window) of the values of the x component of the electric field at that receiver position. +* ``Ey`` is an array containing the time history (for the model time window) of the values of the y component of the electric field at that receiver position. +* ``Ez`` is an array containing the time history (for the model time window) of the values of the z component of the electric field at that receiver position. +* ``Hx`` is an array containing the time history (for the model time window) of the values of the x component of the magnetic field at that receiver position. +* ``Hy`` is an array containing the time history (for the model time window) of the values of the y component of the magnetic field at that receiver position. +* ``Hz`` is an array containing the time history (for the model time window) of the values of the z component of the magnetic field at that receiver position. + +Within each individual ``tx`` group is the following dataset: + +* ``Position`` is the x, y, z position (in metres) of the receiver in the model. + + +Viewing output +============== + +There are a number of free tools available to read HDF5 files. Also MATLAB has high- and low-level functions for reading and writing HDF5 files, i.e. ``h5info`` and ``h5disp`` are useful for returning information and displaying the contents of HDF5 files respectively. gprMax includes some Python modules (in the ``tools`` package) to help you view output data: + +A-scans +------- + +* Plot A-scans using the Python module ``plot_hdf5_Ascan.py``. The module uses matplotlib to plot the time history for the electric and magnetic field components for all receivers in a model (each receiver gets a separate figure window). Usage (from the top-level gprMax directory) is: ``python -m tools.plot_hdf5_Ascan my_outputfile.out``. + +* Plot A-scans using the MATLAB script ``plot_hdf5_Ascan.m``. The script plots the time history for the electric and magnetic field components for all receivers in a model (each receiver gets a separate figure window). + +B-scans +------- + +gprMax produces a separate output file for each trace (A-scan) in the B-scan. + +* Combine the separate output files into one file using the Python module ``outputfiles_mergehdf5.py``. Usage (from the top-level gprMax directory) is: ``python -m tools.outputfiles_mergehdf5 basefilename modelruns``, where ``basefilename`` is the base name file of the output file series, e.g. for ``myoutput1.out``, ``myoutput2.out`` the base file name would be ``myoutput``, and ``modelruns`` is the number of output files to combine. +* Plot an image of the B-scan using the Python module ``plot_hdf5_Bscan.py``. Usage (from the top-level gprMax directory) is: ``python -m tools.plot_hdf5_Bscan my_outputfile.out field``, where ``field`` is the name of field to plot, e.g. ``Ex``, ``Ey`` or ``Ez``. + +