diff --git a/docs/source/index.rst b/docs/source/index.rst index 185c3986..abb868e6 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -29,6 +29,7 @@ gprMax User Guide :maxdepth: 2 :caption: Advanced topics + api python_scripting openmp_mpi gpu diff --git a/docs/source/python_scripting.rst b/docs/source/python_scripting.rst index 3c782788..e11e959e 100644 --- a/docs/source/python_scripting.rst +++ b/docs/source/python_scripting.rst @@ -4,6 +4,8 @@ Scripting the input file ************************ +INTEGRATE THIS WHICH API DOCS + 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. Constants/variables @@ -21,20 +23,3 @@ You can access the following built-in variables from your Python code: * ``current_model_run`` is the current run number of the model that is been executed. * ``inputfile`` is the path and name of the input file. * ``number_model_runs`` is the total number of runs specified when the model was initially executed, i.e. from ``python -m gprMax my_input_file -n number_of_model_runs`` - - -Functions for input commands -============================ - -To make it easier to create commands within a block of Python code, there is a built-in module which contains some of the most commonly used input commands in functional form. For example, to use Python to generate a series of cylinders in a model: - -.. code-block:: none - - #python: - from gprMax.input_cmd_funcs import * - domain = domain(0.2, 0.2, 0.2) - for x in range(0, 8): - cylinder(0.02 + x * 0.02, 0.05, 0, 0.020 + x * 0.02, 0.05, domain[2], 0.005, 'pec’) - #end_python: - -The ``domain`` function will print the ``#domain`` command to the input file and return a variable with the extent of the domain that can be used elsewhere in a Python code block, e.g. in this case with the ``cylinder`` function. The ``cylinder`` function is just a functional version of the ``#cylinder`` command which prints it to the input file. diff --git a/gprMax/__init__.py b/gprMax/__init__.py index d5f536c4..92526be3 100644 --- a/gprMax/__init__.py +++ b/gprMax/__init__.py @@ -7,14 +7,6 @@ Electromagnetic wave propagation simulation software. """ -import os -import sys -# There is a bug with threading and MKL on macOS -# (https://github.com/gprMax/gprMax/issues/195) . Setting the MKL threading -# layer to sequential solves it, but must be done before numpy is imported. -#if sys.platform == 'darwin': -# os.environ["MKL_THREADING_LAYER"] = 'sequential' - from ._version import __version__ from .cmds_single_use import Discretisation from .cmds_single_use import Domain @@ -46,7 +38,6 @@ from .cmds_multiple import PMLCFS from .subgrids.user_objects import SubGridHSG from .subgrids.user_objects import ReferenceRx -# import geometry commands from .cmds_geometry.edge import Edge from .cmds_geometry.plate import Plate from .cmds_geometry.triangle import Triangle diff --git a/gprMax/_version.py b/gprMax/_version.py index e261f36e..1a33f626 100644 --- a/gprMax/_version.py +++ b/gprMax/_version.py @@ -1,4 +1,4 @@ # This is where the version number is set and read by setup.py and conf.py (for the docs) -__version__ = '3.1.6b0' -codename = 'Big Smoke' +__version__ = '3.2.0b0' +codename = '' diff --git a/gprMax/cmds_geometry/add_grass.py b/gprMax/cmds_geometry/add_grass.py index 88c8fd3c..cdaab566 100644 --- a/gprMax/cmds_geometry/add_grass.py +++ b/gprMax/cmds_geometry/add_grass.py @@ -1,14 +1,31 @@ -"""Class for add_grass command.""" +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import numpy as np +from tqdm import tqdm + +import gprMax.config as config from .cmds_geometry import UserObjectGeometry from ..exceptions import CmdInputError -from ..utilities import round_value -from ..materials import Material from ..fractals import FractalSurface from ..fractals import Grass -import gprMax.config as config - -from tqdm import tqdm -import numpy as np +from ..materials import Material +from ..utilities import round_value class AddGrass(UserObjectGeometry): diff --git a/gprMax/cmds_geometry/add_surface_roughness.py b/gprMax/cmds_geometry/add_surface_roughness.py index db5c5ed0..0e8d815e 100644 --- a/gprMax/cmds_geometry/add_surface_roughness.py +++ b/gprMax/cmds_geometry/add_surface_roughness.py @@ -1,12 +1,30 @@ -"""Class for surface roughness command.""" +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import numpy as np +from tqdm import tqdm + +import gprMax.config as config from .cmds_geometry import UserObjectGeometry +from ..cython.geometry_primitives import build_box from ..exceptions import CmdInputError from ..fractals import FractalSurface from ..utilities import round_value -import gprMax.config as config - -from tqdm import tqdm -import numpy as np class AddSurfaceRoughness(UserObjectGeometry): diff --git a/gprMax/cmds_geometry/add_surface_water.py b/gprMax/cmds_geometry/add_surface_water.py index 4c381338..654e566c 100644 --- a/gprMax/cmds_geometry/add_surface_water.py +++ b/gprMax/cmds_geometry/add_surface_water.py @@ -1,12 +1,29 @@ -"""Class for surface water command.""" -from .cmds_geometry import UserObjectGeometry -from ..exceptions import CmdInputError -from ..utilities import round_value -from ..materials import Material -import gprMax.config as config +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . from tqdm import tqdm +import gprMax.config as config +from .cmds_geometry import UserObjectGeometry +from ..exceptions import CmdInputError +from ..materials import Material +from ..utilities import round_value + class AddSurfaceWater(UserObjectGeometry): """Allows you to add surface water to a :class:`gprMax.cmds_geometry.fractal_box.FractalBox` in the model. diff --git a/gprMax/cmds_geometry/box.py b/gprMax/cmds_geometry/box.py index 88db127f..8020fc6e 100644 --- a/gprMax/cmds_geometry/box.py +++ b/gprMax/cmds_geometry/box.py @@ -1,12 +1,29 @@ -"""Class for triangle command.""" +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import numpy as np +from tqdm import tqdm + +import gprMax.config as config from .cmds_geometry import UserObjectGeometry +from ..cython.geometry_primitives import build_box from ..exceptions import CmdInputError from ..materials import Material -from ..cython.geometry_primitives import build_box -import gprMax.config as config - -from tqdm import tqdm -import numpy as np class Box(UserObjectGeometry): diff --git a/gprMax/cmds_geometry/build_templates.py b/gprMax/cmds_geometry/build_templates.py index 54ea0547..a356aee7 100644 --- a/gprMax/cmds_geometry/build_templates.py +++ b/gprMax/cmds_geometry/build_templates.py @@ -1,4 +1,25 @@ -from jinja2 import Environment, PackageLoader, select_autoescape +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +from jinja2 import Environment +from jinja2 import PackageLoader +from jinja2 import select_autoescape + env = Environment( loader=PackageLoader(__name__, 'templates'), ) diff --git a/gprMax/cmds_geometry/cmds_geometry.py b/gprMax/cmds_geometry/cmds_geometry.py index 2f56fc2f..5847abb2 100644 --- a/gprMax/cmds_geometry/cmds_geometry.py +++ b/gprMax/cmds_geometry/cmds_geometry.py @@ -1,3 +1,21 @@ +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + class UserObjectGeometry: """Specific Geometry object.""" diff --git a/gprMax/cmds_geometry/cylinder.py b/gprMax/cmds_geometry/cylinder.py index 30d2fcb3..0c43d6fe 100644 --- a/gprMax/cmds_geometry/cylinder.py +++ b/gprMax/cmds_geometry/cylinder.py @@ -1,12 +1,29 @@ -"""Class for cylinder command.""" +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import numpy as np +from tqdm import tqdm + +import gprMax.config as config from .cmds_geometry import UserObjectGeometry +from ..cython.geometry_primitives import build_cylinder from ..exceptions import CmdInputError from ..materials import Material -from ..cython.geometry_primitives import build_cylinder - -from tqdm import tqdm -import numpy as np -import gprMax.config as config class Cylinder(UserObjectGeometry): diff --git a/gprMax/cmds_geometry/cylindrical_sector.py b/gprMax/cmds_geometry/cylindrical_sector.py index 8cc5e8e4..05781e59 100644 --- a/gprMax/cmds_geometry/cylindrical_sector.py +++ b/gprMax/cmds_geometry/cylindrical_sector.py @@ -1,13 +1,29 @@ -"""Class for cylinder command.""" +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import numpy as np +from tqdm import tqdm + +import gprMax.config as config from .cmds_geometry import UserObjectGeometry +from ..cython.geometry_primitives import build_cylindrical_sector from ..exceptions import CmdInputError from ..materials import Material -from ..cython.geometry_primitives import build_cylindrical_sector - -from tqdm import tqdm -import numpy as np -import gprMax.config as config - class CylindricalSector(UserObjectGeometry): diff --git a/gprMax/cmds_geometry/edge.py b/gprMax/cmds_geometry/edge.py index d5950755..3d01eb01 100644 --- a/gprMax/cmds_geometry/edge.py +++ b/gprMax/cmds_geometry/edge.py @@ -1,12 +1,29 @@ -"""Class for edge command.""" +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +from tqdm import tqdm + +import gprMax.config as config from .cmds_geometry import UserObjectGeometry -from ..exceptions import CmdInputError from ..cython.geometry_primitives import build_edge_x from ..cython.geometry_primitives import build_edge_y from ..cython.geometry_primitives import build_edge_z -import gprMax.config as config - -from tqdm import tqdm +from ..exceptions import CmdInputError class Edge(UserObjectGeometry): diff --git a/gprMax/cmds_geometry/fractal_box.py b/gprMax/cmds_geometry/fractal_box.py index 3d346564..913f10c5 100644 --- a/gprMax/cmds_geometry/fractal_box.py +++ b/gprMax/cmds_geometry/fractal_box.py @@ -1,11 +1,28 @@ -"""Class for surface roughness command.""" +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import numpy as np +from tqdm import tqdm + +import gprMax.config as config from .cmds_geometry import UserObjectGeometry from ..exceptions import CmdInputError from ..fractals import FractalVolume -import gprMax.config as config - -from tqdm import tqdm -import numpy as np class FractalBox(UserObjectGeometry): diff --git a/gprMax/cmds_geometry/fractal_box_builder.py b/gprMax/cmds_geometry/fractal_box_builder.py index b61c404b..73313f1e 100644 --- a/gprMax/cmds_geometry/fractal_box_builder.py +++ b/gprMax/cmds_geometry/fractal_box_builder.py @@ -1,15 +1,32 @@ -"""Class for fractal box modifications command.""" -from .cmds_geometry import UserObjectGeometry -from ..exceptions import CmdInputError -from ..config import dtypes -from ..cython.geometry_primitives import build_voxels_from_array -from ..cython.geometry_primitives import build_voxels_from_array_mask +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . import numpy as np +from .cmds_geometry import UserObjectGeometry +from ..config import dtypes +from ..cython.geometry_primitives import build_voxels_from_array +from ..cython.geometry_primitives import build_voxels_from_array_mask +from ..exceptions import CmdInputError floattype = dtypes['float_or_double'] + class FractalBoxBuilder(UserObjectGeometry): """Internal class for fractal box modifications. This class should be used internally only when surface modification have been made to a fractal box""" diff --git a/gprMax/cmds_geometry/geometry_objects_read.py b/gprMax/cmds_geometry/geometry_objects_read.py index a32dd161..034d10e0 100644 --- a/gprMax/cmds_geometry/geometry_objects_read.py +++ b/gprMax/cmds_geometry/geometry_objects_read.py @@ -1,8 +1,28 @@ +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import os + +from tqdm import tqdm + +import gprMax.config as config from .cmds_geometry import UserObjectGeometry from ..exceptions import CmdInputError -import os -from tqdm import tqdm -import gprMax.config as config class GeometryObjectsRead(UserObjectGeometry): diff --git a/gprMax/cmds_geometry/plate.py b/gprMax/cmds_geometry/plate.py index a68dabfe..9852e02c 100644 --- a/gprMax/cmds_geometry/plate.py +++ b/gprMax/cmds_geometry/plate.py @@ -1,12 +1,29 @@ -"""Class for edge command.""" +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +from tqdm import tqdm + +import gprMax.config as config from .cmds_geometry import UserObjectGeometry -from ..exceptions import CmdInputError from ..cython.geometry_primitives import build_face_yz from ..cython.geometry_primitives import build_face_xz from ..cython.geometry_primitives import build_face_xy -import gprMax.config as config - -from tqdm import tqdm +from ..exceptions import CmdInputError class Plate(UserObjectGeometry): diff --git a/gprMax/cmds_geometry/sphere.py b/gprMax/cmds_geometry/sphere.py index 1d568744..80e806cd 100644 --- a/gprMax/cmds_geometry/sphere.py +++ b/gprMax/cmds_geometry/sphere.py @@ -1,12 +1,29 @@ -"""Class for sphere command.""" +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import numpy as np +from tqdm import tqdm + +import gprMax.config as config from .cmds_geometry import UserObjectGeometry +from ..cython.geometry_primitives import build_sphere from ..exceptions import CmdInputError from ..materials import Material -from ..cython.geometry_primitives import build_sphere -import gprMax.config as config - -from tqdm import tqdm -import numpy as np class Sphere(UserObjectGeometry): diff --git a/gprMax/cmds_geometry/triangle.py b/gprMax/cmds_geometry/triangle.py index ecac6380..37086278 100644 --- a/gprMax/cmds_geometry/triangle.py +++ b/gprMax/cmds_geometry/triangle.py @@ -1,12 +1,29 @@ -"""Class for triangle command.""" +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +import numpy as np +from tqdm import tqdm + +import gprMax.config as config from .cmds_geometry import UserObjectGeometry +from ..cython.geometry_primitives import build_triangle from ..exceptions import CmdInputError from ..materials import Material -from ..cython.geometry_primitives import build_triangle -import gprMax.config as config - -from tqdm import tqdm -import numpy as np class Triangle(UserObjectGeometry): diff --git a/gprMax/cmds_multiple.py b/gprMax/cmds_multiple.py index 98852f44..f67ff313 100644 --- a/gprMax/cmds_multiple.py +++ b/gprMax/cmds_multiple.py @@ -15,38 +15,36 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . +import sys +import numpy as np +from tqdm import tqdm -"""Object which can be created multiple times.""" import gprMax.config as config from .config import z0 from .config import dtypes -from .utilities import round_value + from .cmds_geometry.cmds_geometry import UserObjectGeometry -from .waveforms import Waveform as WaveformUser +from .exceptions import CmdInputError +from .geometry_outputs import GeometryObjects as GeometryObjectsUser +from .materials import Material as MaterialUser +from .materials import PeplinskiSoil as PeplinskiSoilUser +from .pml import CFSParameter +from .pml import CFS +from .receivers import Rx as RxUser +from .snapshots import Snapshot as SnapshotUser from .sources import VoltageSource as VoltageSourceUser from .sources import HertzianDipole as HertzianDipoleUser from .sources import MagneticDipole as MagneticDipoleUser from .sources import TransmissionLine as TransmissionLineUser -from .snapshots import Snapshot as SnapshotUser -from .receivers import Rx as RxUser -from .materials import Material as MaterialUser -from .materials import PeplinskiSoil as PeplinskiSoilUser -from .geometry_outputs import GeometryObjects as GeometryObjectsUser -from .pml import CFSParameter -from .pml import CFS from .subgrids.base import SubGridBase - -from .exceptions import CmdInputError - - -import numpy as np -from tqdm import tqdm +from .utilities import round_value +from .waveforms import Waveform as WaveformUser floattype = dtypes['float_or_double'] class UserObjectMulti: - """Specific multiobject object.""" + """Object that can occur multiple times in a model.""" def __init__(self, **kwargs): """Constructor.""" @@ -58,7 +56,7 @@ class UserObjectMulti: def __str__(self): - """Readble user string as per hash commands.""" + """Readable user string as per hash commands.""" s = '' for k, v in self.kwargs.items(): if isinstance(v, tuple) or isinstance(v, list): @@ -86,6 +84,7 @@ class Waveform(UserObjectMulti): :param freq: The centre frequency of the waveform (Hertz) :type freq: float, non-optional """ + def __init__(self, **kwargs): super().__init__(**kwargs) self.order = 0 @@ -136,6 +135,7 @@ class VoltageSource(UserObjectMulti): :param stop: Time to remove the source :type stop: float, optional """ + def __init__(self, **kwargs): super().__init__(**kwargs) self.order = 1 @@ -296,7 +296,6 @@ class HertzianDipole(UserObjectMulti): except KeyError: h.start = 0 h.stop = grid.timewindow - print(grid.timewindow) startstop = ' ' h.calculate_waveform_values(grid) @@ -414,6 +413,7 @@ class TransmissionLine(UserObjectMulti): :param stop: Time to remove the source :type stop: float, optional """ + def __init__(self, **kwargs): """Constructor.""" super().__init__(**kwargs) @@ -505,7 +505,6 @@ class Rx(UserObjectMulti): :param outputs: is a list of outputs with this receiver. It can be any selection from Ex, Ey, Ez, Hx, Hy, Hz, Ix, Iy, or Iz. :type outputs: list, non-optional - """ def __init__(self, **kwargs): @@ -554,6 +553,7 @@ class Rx(UserObjectMulti): return r + class RxArray(UserObjectMulti): """Provides a simple method of defining multiple output points in the model. @@ -699,6 +699,7 @@ class Snapshot(UserObjectMulti): grid.snapshots.append(s) + class Material(UserObjectMulti): """Allows you to introduce a material into the model described by a set of constitutive parameters. @@ -1219,6 +1220,7 @@ class PMLCFS(UserObjectMulti): grid.cfs.append(cfs) + class Subgrid(UserObjectMulti): def __init__(self, **kwargs): @@ -1234,6 +1236,7 @@ class Subgrid(UserObjectMulti): else: raise Exception('This Object is Unknown to gprMax') + class SubgridHSG(UserObjectMulti): def __init__(self, **kwargs): diff --git a/gprMax/cmds_single_use.py b/gprMax/cmds_single_use.py index 62d3f2a6..82de288d 100644 --- a/gprMax/cmds_single_use.py +++ b/gprMax/cmds_single_use.py @@ -16,10 +16,17 @@ # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . +import decimal as d import inspect import os import sys -import decimal as d + +from colorama import init +from colorama import Fore +from colorama import Style +init() +import numpy as np +from scipy import interpolate import gprMax.config as config from .config import c @@ -29,13 +36,6 @@ from .exceptions import CmdInputError from .waveforms import Waveform from .utilities import round_value -import numpy as np -from scipy import interpolate -from colorama import init -from colorama import Fore -from colorama import Style -init() - floattype = dtypes['float_or_double'] @@ -44,18 +44,17 @@ class Properties: class UserObjectSingle: + """Object that can only occur a single time in a model.""" def __init__(self, **kwargs): - # each single command has an order to specify the order in which - # the commands are constructed. IE. discretisation must be + # Each single command has an order to specify the order in which + # the commands are constructed, e.g. discretisation must be # created before the domain self.order = None self.kwargs = kwargs self.props = Properties() - # auto translate self.autotranslate = True - for k, v in kwargs.items(): setattr(self.props, k, v) @@ -74,8 +73,6 @@ class Domain(UserObjectSingle): """ def __init__(self, **kwargs): - # dont need to define parameters in advance. Just catch errors - # when they occur super().__init__(**kwargs) self.order = 2 @@ -90,8 +87,6 @@ class Domain(UserObjectSingle): return s def create(self, G, uip): - - # code to create the gprMax domain as per input_cmds_singleuse.py try: G.nx, G.ny, G.nz = uip.discretise_point(self.kwargs['p1']) except KeyError: @@ -122,7 +117,6 @@ class Domain(UserObjectSingle): G.dt = 1 / (c * np.sqrt((1 / G.dx) * (1 / G.dx) + (1 / G.dy) * (1 / G.dy) + (1 / G.dz) * (1 / G.dz))) G.mode = '3D' - # Round down time step to nearest float with precision one less than hardware maximum. Avoids inadvertently exceeding the CFL due to binary representation of floating point number. # Round down time step to nearest float with precision one less than hardware maximum. # Avoids inadvertently exceeding the CFL due to binary representation of floating point number. G.dt = round_value(G.dt, decimalplaces=d.getcontext().prec - 1) @@ -164,9 +158,8 @@ class Discretisation(UserObjectSingle): :param p1: Specify discretisation in x, y, z direction :type p1: list of floats, non-optional """ + def __init__(self, **kwargs): - # dont need to define parameters in advance. Just catch errors - # when they occur super().__init__(**kwargs) self.order = 1 @@ -181,11 +174,9 @@ class Discretisation(UserObjectSingle): return s def create(self, G, uip): - try: G.dl = np.array(self.kwargs['p1']) G.dx, G.dy, G.dz = self.kwargs['p1'] - except KeyError: raise CmdInputError('Discretisation requires a point') @@ -208,9 +199,8 @@ class TimeWindow(UserObjectSingle): :param iterations: Required number of iterations :type iterations: int, optional """ + def __init__(self, **kwargs): - # dont need to define parameters in advance. Just catch errors - # when they occur super().__init__(**kwargs) self.order = 4 @@ -263,8 +253,6 @@ class Messages(UserObjectSingle): """ def __init__(self, **kwargs): - # dont need to define parameters in advance. Just catch errors - # when they occur super().__init__(**kwargs) self.order = 0 @@ -296,8 +284,6 @@ class Title(UserObjectSingle): """ def __init__(self, **kwargs): - # dont need to define parameters in advance. Just catch errors - # when they occur super().__init__(**kwargs) self.order = 5 @@ -312,6 +298,7 @@ class Title(UserObjectSingle): if config.is_messages(): print('Model title: {}'.format(G.title)) + class NumThreads(UserObjectSingle): """Allows you to control how many OpenMP threads (usually the number of physical CPU cores available) are used when running the model. @@ -319,9 +306,8 @@ class NumThreads(UserObjectSingle): :param n: Number of threads. :type n: int, optional """ + def __init__(self, **kwargs): - # dont need to define parameters in advance. Just catch errors - # when they occur super().__init__(**kwargs) self.order = 6 @@ -356,28 +342,24 @@ class NumThreads(UserObjectSingle): print('GPU solving using: {} - {}'.format(G.gpu.deviceID, G.gpu.name)) -# Time step stability factor class TimeStepStabilityFactor(UserObjectSingle): """Factor by which to reduce the time step from the CFL limit. :param f: Factor to multiple time step. :type f: float, optional """ + def __init__(self, **kwargs): - # dont need to define parameters in advance. Just catch errors - # when they occur super().__init__(**kwargs) self.order = 7 def __str__(self): - try: return '#time_step_stability_factor: {}'.format(self.kwargs['f']) except KeyError: return '#time_step_stability_factor:' def create(self, G, uip): - try: f = self.kwargs['f'] except KeyError: @@ -410,14 +392,12 @@ class PMLCells(UserObjectSingle): :param zmax: Thickness of PML on top side. :type zmax: int, optional """ + def __init__(self, **kwargs): - # dont need to define parameters in advance. Just catch errors - # when they occur super().__init__(**kwargs) self.order = 8 def create(self, G, uip): - try: thickness = self.kwargs['thickness'] @@ -452,8 +432,6 @@ class SrcSteps(UserObjectSingle): """ def __init__(self, **kwargs): - # dont need to define parameters in advance. Just catch errors - # when they occur super().__init__(**kwargs) self.order = 9 @@ -462,7 +440,7 @@ class SrcSteps(UserObjectSingle): G.srcsteps = uip.discretise_point(self.kwargs['p1']) except KeyError: raise CmdInputError('#src_steps: requires exactly three parameters') - # src_steps + if config.is_messages(): print('Simple sources will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.srcsteps[0] * G.dx, G.srcsteps[1] * G.dy, G.srcsteps[2] * G.dz)) @@ -475,8 +453,6 @@ class RxSteps(UserObjectSingle): """ def __init__(self, **kwargs): - # dont need to define parameters in advance. Just catch errors - # when they occur super().__init__(**kwargs) self.order = 10 @@ -485,6 +461,7 @@ class RxSteps(UserObjectSingle): G.rxsteps = uip.discretise_point(self.kwargs['p1']) except KeyError: raise CmdInputError('#rx_steps: requires exactly three parameters') + if config.is_messages(): print('All receivers will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.rxsteps[0] * G.dx, G.rxsteps[1] * G.dy, G.rxsteps[2] * G.dz)) @@ -502,7 +479,6 @@ class ExcitationFile(UserObjectSingle): """ def create(self, G, uip): - # Excitation file for user-defined source waveforms try: kwargs = dict() excitationfile = self.kwargs['filepath'] @@ -574,8 +550,6 @@ class OutputDir(UserObjectSingle): :type dir: str, non-optional """ def __init__(self, **kwargs): - # dont need to define parameters in advance. Just catch errors - # when they occur super().__init__(**kwargs) self.order = 11 diff --git a/gprMax/config.py b/gprMax/config.py index 7328cd74..7bce848b 100644 --- a/gprMax/config.py +++ b/gprMax/config.py @@ -17,22 +17,21 @@ # along with gprMax. If not, see . import os - -import cython -import numpy as np -from scipy.constants import c -from scipy.constants import mu_0 as m0 -from scipy.constants import epsilon_0 as e0 -import sys - -from .utilities import get_terminal_width -from .utilities import get_host_info from pathlib import Path +import sys from colorama import init from colorama import Fore from colorama import Style init() +import cython +import numpy as np +from scipy.constants import c +from scipy.constants import mu_0 as m0 +from scipy.constants import epsilon_0 as e0 + +from .utilities import get_host_info +from .utilities import get_terminal_width # Impedance of free space (Ohms) @@ -51,12 +50,10 @@ z0 = np.sqrt(m0 / e0) general = {'inputfilepath': 'gprMax', 'outputfilepath': 'gprMax', 'messages': True, 'progressbars': True, 'mode': '3D', 'cpu': True, 'cuda': False, 'opencl': False, 'autotranslate': False} - def is_messages(): """Function to return messages.""" return general['messages'] - # Store information about host machine hostinfo = get_host_info() @@ -98,38 +95,39 @@ elif precision == 'double': class ModelConfig(): + """Configuration parameters for a model. + N.B. Multiple models can exist within a simulation + """ def __init__(self, sim_config, i): + """ + Args: + sim_config (SimConfig): Simulation level configuration object. + i (int): Current model number. + """ + self.sim_config = sim_config self.reuse_geometry = False - - # current model number (indexed from 0) - self.i = i - - parts = self.sim_config.output_file_path.parts + self.i = i # Indexed from 0 if not sim_config.single_model: - # 1 indexed - self.appendmodelnumber = str(self.i + 1) + self.appendmodelnumber = str(self.i + 1) # Indexed from 1 else: self.appendmodelnumber = '' - # outputfilepath for specific model - self.output_file_path = Path(*parts[:-2], parts[-1] + self.appendmodelnumber).with_suffix('') + # Output file path for specific model + parts = self.sim_config.output_file_path.with_suffix('').parts + self.output_file_path = Path(*parts[:-1], parts[-1] + self.appendmodelnumber) self.output_file_path_ext = self.output_file_path.with_suffix('.out') - # make a snapshot directory + # Make a snapshot directory self.snapshot_dir = '_snaps' + # String to print at start of each model run inputfilestr_f = '\n--- Model {}/{}, input file: {}' self.inputfilestr = inputfilestr_f.format(self.i + 1, self.sim_config.model_end, self.sim_config.input_file_path) - # string to print at start of each model run self.next_model = Fore.GREEN + '{} {}\n'.format(self.inputfilestr, '-' * (get_terminal_width() - 1 - len(self.inputfilestr))) + Style.RESET_ALL - # Add the current model run to namespace that can be accessed by - # user in any Python code blocks in input file - #self.usernamespace['current_model_run'] = self.i + 1 - def get_scene(self): if self.sim_config.scenes: return self.sim_config.scenes[self.i] @@ -146,24 +144,15 @@ class ModelConfig(): class SimulationConfig: + """Configuration parameters for a standard simulation. + N.B. A simulation can consist of multiple models. + """ def __init__(self, args): - """Adapter for args into Simulation level configuration""" - - # adapt the arg properties to link smoothly with MPIRunner(), CPURunner() etc.. - - # args.inputfile - # args.n - # args.task - # args.restart - # args.mpi - # args.mpi_no_spawn - # args.mpicomm - # args.gpu - # args.benchmark - # args.geometry_only - # args.geometry_fixed - # args.write_processed + """ + Args: + args (Namespace): Arguments from either API or CLI. + """ self.args = args self.n_models = args.n @@ -177,20 +166,19 @@ class SimulationConfig: self.geometry_only = args.geometry_only self.write_processed = args.write_processed - # subgrid parameter may not exist if user uses CLI api + # Subgrid parameter may not exist if user enters via CLI try: self.subgrid = args.subgrid except AttributeError: - # this must be CLI user. No subgrids are available self.subgrid = False - # scenes parameter may not exist if user uses CLI api + # Scenes parameter may not exist if user enters via CLI try: self.scenes = args.scenes except AttributeError: self.scenes = [] - # set more complex parameters + # Set more complex parameters self.set_input_file_path() self.set_output_file_path() self.set_model_start_end() @@ -202,10 +190,8 @@ class SimulationConfig: else: self.single_model = False - # for example def set_model_start_end(self): - - # set range for number of models to run (internally 0 index) + # Set range for number of models to run (internally 0 index) if self.args.task: # Job array feeds args.n number of single tasks modelstart = self.args.task - 1 @@ -224,17 +210,15 @@ class SimulationConfig: pass def set_input_file_path(self): - """Function to set to inputfile path""" - - # if the API is in use an id for the simulation must be provided. + # If the API is in use an id for the simulation must be provided. if self.args.inputfile is None: self.input_file_path = Path(self.args.outputfile) else: self.input_file_path = Path(self.args.inputfile) def set_output_file_path(self): - # output file can be provided by the user. if they havent provided None - # use the inputfilefile path instead + # Output file path can be provided by the user. If they havent provided one + # use the inputfile file path instead try: self.output_file_path = Path(self.args.outputfile) except AttributeError: @@ -242,6 +226,9 @@ class SimulationConfig: class SimulationConfigMPI(SimulationConfig): + """Configuration parameters for a MPI simulation. + N.B. A simulation can consist of multiple models. + """ def __init__(self, args): super().__init__(args) @@ -253,14 +240,33 @@ class SimulationConfigMPI(SimulationConfig): def create_simulation_config(args): + """Create simulation level configuration object to hold simulation + level parameters. + + Args: + args (Namespace): Arguments from either API or CLI. + + Returns: + sc (SimulationConfig): Simulation configuration object. + """ if not args.mpi and not args.mpi_no_spawn: sc = SimulationConfig(args) elif args.mpi: sc = SimulationConfigMPI(args) + return sc def create_model_config(sim_config, i): + """Create model level configuration object to hold model level + parameters. + + Args: + sim_config (SimConfig): Simulation level configuration object. + i (int): Current model number. + """ + mc = ModelConfig(sim_config, i) + return mc diff --git a/gprMax/contexts.py b/gprMax/contexts.py index 04842070..ed2e7f1e 100644 --- a/gprMax/contexts.py +++ b/gprMax/contexts.py @@ -15,28 +15,30 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -from .utilities import get_terminal_width -from .utilities import timer -from .model_build_run import ModelBuildRun + import datetime + from .config import create_model_config +from .model_build_run import ModelBuildRun from .solvers import create_solver from .solvers import create_G +from .utilities import get_terminal_width +from .utilities import timer class Context(): + """Generic context for the model to run in. Sub-class with specific contexts + e.g. an MPI context. + """ def __init__(self, sim_config): - """Context for the model to run in. Sub-class this with contexts - i.e. an MPI context. - + """ Args: sim_config (SimConfig): Simulation level configuration object. - solver (Solver): FDTD general solver object. """ + self.sim_config = sim_config - self.model_range = range(sim_config.model_start, - sim_config.model_end) + self.model_range = range(sim_config.model_start, sim_config.model_end) self.tsimend = 0 self.tsimstart = 1 @@ -52,11 +54,14 @@ class Context(): print(s) def make_time_report(self): - """Function to generate a string for the total simulation time bas""" + """Function to generate a string for the total simulation time.""" pass class NoMPIContext(Context): + """Standard context - models are run one after another and each model + is parallelised using either OpenMP (CPU) or CUDA (GPU). + """ def _run(self): """Specialise how the models are farmed out.""" @@ -64,9 +69,9 @@ class NoMPIContext(Context): for i in self.model_range: model_config = create_model_config(self.sim_config, i) - # always create a solver for the first model - # the next model to run only gets a new solver if the - # geometry is not re used. + # Always create a solver for the first model. + # The next model to run only gets a new solver if the + # geometry is not re-used. if i != 0 and self.sim_config.geometry_fixed: model_config.reuse_geometry = True else: @@ -82,7 +87,9 @@ class NoMPIContext(Context): def make_time_report(self): """Function to specialise the time reporting for the standard Simulation - context.""" + context. + """ + sim_time = datetime.timedelta(seconds=self.tsimend - self.tsimstart) s = '\n=== Simulation on {} completed in [HH:MM:SS]: {}' s = s.format(self.simconfig.hostinfo['hostname'], sim_time) @@ -90,6 +97,15 @@ class NoMPIContext(Context): class MPIContext(Context): + """Mixed mode MPI/OpenMP/CUDA context - MPI task farm is used to distribute + models, and each model parallelised using either OpenMP (CPU) + or CUDA (GPU). + """ + + def __init__(self): + super().__init__() + + from mpi4py import MPI def _run(self): pass @@ -108,10 +124,18 @@ class MPINoSpawnContext(Context): def create_context(sim_config): - """Create a context in which to run the simulation. i.e MPI.""" - if sim_config.mpi_no_spawn: + """Create a context in which to run the simulation. i.e MPI. + + Args: + sim_config (SimConfig): Simulation level configuration object. + + Returns: + context (Context): Context for the model to run in. + """ + + if sim_config.mpi: context = MPIContext(sim_config) - elif sim_config.mpi: + elif sim_config.mpi_no_spawn: context = MPINoSpawnContext(sim_config) else: context = NoMPIContext(sim_config) diff --git a/gprMax/fields_outputs.py b/gprMax/fields_outputs.py index 26a5212f..decaab96 100644 --- a/gprMax/fields_outputs.py +++ b/gprMax/fields_outputs.py @@ -31,14 +31,12 @@ def store_outputs(G): iteration (int): Current iteration number. Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic field values. - G (class): Grid class instance - holds essential parameters describing - the model. + G (Grid): Holds essential parameters describing the model. """ iteration = G.iteration - #import pdb; pdb.set_trace() - Ex, Ey, Ez, Hx, Hy, Hz = G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz + for rx in G.rxs: for output in rx.outputs: # Store electric or magnetic field components @@ -107,8 +105,7 @@ def write_hdf5_main_grid_outputfile(outputfile, G): Args: outputfile (str): Name of the output file. - Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic field values. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ write_data(outputfile, G) @@ -119,8 +116,7 @@ def write_hdf5_sub_grid_outputfile(outputfile, G): Args: outputfile (str): Name of the output file. - Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic field values. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ stem = outputfile.stem @@ -129,13 +125,13 @@ def write_hdf5_sub_grid_outputfile(outputfile, G): for sg in G.subgrids: - # create an outputfile for each subgrid + # Create an outputfile for each subgrid fp = stem + '_' + sg.name + suffix fp = parent / Path(fp) f = write_data(fp, sg) - # write some additional meta data about the subgrid + # Write some additional meta data about the subgrid f.attrs['is_os_sep'] = sg.is_os_sep f.attrs['pml_separation'] = sg.pml_separation f.attrs['subgrid_pml_thickness'] = sg.pmlthickness['x0'] @@ -149,10 +145,7 @@ def write_data(outputfile, G): Args: outputfile (str): Name of the output file. - Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic - field values. - G (class): Grid class instance - holds essential parameters describing - the model. + G (Grid): Holds essential parameters describing the model. """ f = h5py.File(outputfile, 'w') diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index b297fe00..fcbc6375 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.py @@ -16,11 +16,11 @@ # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . +import argparse + from .config import create_simulation_config from .contexts import create_context -from .solvers import create_solver -import argparse def run( scenes=None, @@ -34,41 +34,55 @@ def run( mpicomm=None, gpu=None, subgrid=None, - benchmark=False, geometry_only=False, geometry_fixed=False, write_processed=False, ): - """Run the simulation for the given list of scenes. + """This is the main function for gprMax when entering as application + programming interface (API). Run the simulation for the + given list of scenes. - :param scenes: List of the scenes to run the model. Multiple scene objects can given in order to run multiple simulation runs. Each scene must contain the essential simulation objects + :param scenes: list of the scenes to run the model. Multiple scene objects can given in order to run multiple simulation runs. Each scene must contain the essential simulation objects :type scenes: list, optional - :param inputfile: Input file path. Can also run simulation by providing an input file. + + :param inputfile: input file path. Can also run simulation by providing an input file. :type inputfile: str, optional - :param outputfile: File path to the output data file. + + :param outputfile: file path to the output data file. :type outputfile: str, non-optional - :param n: Number of required simulation runs. + + :param n: number of required simulation runs. :type n: int, non-optional + :param task: task identifier (model number) when running simulation as a job array on open grid scheduler/grid engine. for further details see the parallel performance section of the user guide :type task: int, optional + :param restart: model number to start/restart simulation from. It would typically be used to restart a series of models from a specific model number, with the n argument, e.g. to restart from A-scan 45 when creating a B-scan with 60 traces :type restart: int, optional + :param mpi: number of Message Passing Interface (MPI) tasks, i.e. master + workers, for MPI task farm. This option is most usefully combined with n to allow individual models to be farmed out using a MPI task farm, e.g. to create a B-scan with 60 traces and use MPI to farm out each trace1. For further details see the parallel performance section of the User Guide :type mpi: int, optional - :param mpi_no_spawn: use MPI task farm without spawn mechanism. For further details see the parallel performance section of the User Guide. + + :param mpi_no_spawn: flag to use MPI task farm without spawn mechanism. For further details see the parallel performance section of the User Guide. :type mpi_no_spawn: bool, optional - :param gpu: Flag to use NVIDIA GPU or list of NVIDIA GPU device ID(s) for specific GPU card(s) + + :param gpu: flag to use NVIDIA GPU or list of NVIDIA GPU device ID(s) for specific GPU card(s) :type gpu: list or bool, optional - :param subgrid: Use sub-gridding. + + :param subgrid: flag to use sub-gridding. :type subgrid: bool, optional - :param benchmark: Switch on benchmarking mode. This can be used to benchmark the threading (parallel) performance of gprMax on different hardware. For further details see the benchmarking section of the User Guide - :type benchmark: bool, optional + :param geometry_only: build a model and produce any geometry views but do not run the simulation. :type geometry_only: bool, optional + :param geometry_fixed: build a model and produce any geometry views but do not run the simulation. :type geometry_fixed: bool, optional + :param write_processed: write another input file after any Python code and + in the original input file has been processed. + :type write_processed: bool, optional """ + class ImportArguments: pass @@ -85,7 +99,6 @@ def run( args.mpicomm = mpicomm args.gpu = gpu args.subgrid=subgrid - args.benchmark = benchmark args.geometry_only = geometry_only args.geometry_fixed = geometry_fixed args.write_processed = write_processed @@ -94,7 +107,7 @@ def run( def main(): - """This is the main function for gprMax.""" + """Main function for gprMax when entering from the command line interface (CLI).""" # Parse command line arguments parser = argparse.ArgumentParser(prog='gprMax', formatter_class=argparse.ArgumentDefaultsHelpFormatter) @@ -106,7 +119,6 @@ def main(): parser.add_argument('--mpi-no-spawn', action='store_true', default=False, help='flag to use MPI without spawn mechanism') parser.add_argument('--mpi-worker', action='store_true', default=False, help=argparse.SUPPRESS) parser.add_argument('-gpu', type=int, action='append', nargs='*', help='flag to use Nvidia GPU or option to give list of device ID(s)') - parser.add_argument('-benchmark', action='store_true', default=False, help='flag to switch on benchmarking mode') parser.add_argument('--geometry-only', action='store_true', default=False, help='flag to only build model and produce geometry file(s)') parser.add_argument('--geometry-fixed', action='store_true', default=False, help='flag to not reprocess model geometry, e.g. for B-scans where the geometry is fixed') parser.add_argument('--write-processed', action='store_true', default=False, help='flag to write an input file after any Python code and include commands in the original input file have been processed') @@ -114,7 +126,13 @@ def main(): run_main(args) + def run_main(args): + """Called by either run (API) or main (CLI). + + Args: + args (Namespace): arguments from either API or CLI. + """ sim_config = create_simulation_config(args) context = create_context(sim_config) diff --git a/gprMax/grid.py b/gprMax/grid.py index 8dd6ded4..643cf425 100644 --- a/gprMax/grid.py +++ b/gprMax/grid.py @@ -17,6 +17,7 @@ # along with gprMax. If not, see . from collections import OrderedDict +import decimal as d from colorama import init from colorama import Fore @@ -24,6 +25,7 @@ from colorama import Style init() import numpy as np np.seterr(invalid='raise') +from scipy.constants import c import gprMax.config as config from .exceptions import GeneralError @@ -33,10 +35,6 @@ from .utilities import fft_power from .utilities import human_size from .utilities import round_value -import decimal as d -from scipy.constants import c - - class Grid(object): """Generic grid/mesh.""" @@ -104,14 +102,16 @@ class Grid(object): class FDTDGrid(Grid): - """ - Holds attributes associated with the entire grid. A convenient - way for accessing regularly used parameters. + """Holds attributes associated with entire grid. A convenient + way for accessing regularly used parameters. """ def __init__(self): self.title = '' self.memoryusage = 0 + self.name = 'Main' + self.gpu = None + self.outputdirectory = '' self.nx = 0 self.ny = 0 @@ -120,6 +120,7 @@ class FDTDGrid(Grid): self.dy = 0 self.dz = 0 self.dt = 0 + self.iteration = 0 self.iterations = 0 self.timewindow = 0 @@ -148,14 +149,10 @@ class FDTDGrid(Grid): self.rxsteps = [0, 0, 0] self.snapshots = [] self.subgrids = [] - self.gpu = None - self.name = 'Main' - self.outputdirectory = '' - self.iteration = 0 + def initialise_geometry_arrays(self): - """ - Initialise an array for volumetric material IDs (solid); + """Initialise an array for volumetric material IDs (solid); boolean arrays for specifying whether materials can have dielectric smoothing (rigid); and an array for cell edge IDs (ID). Solid and ID arrays are initialised to free_space (one); @@ -286,19 +283,19 @@ class FDTDGrid(Grid): self.ID[2, 0, :, :] = 0 self.ID[2, 1, :, :] = 0 - def tmy(): + def tmy(self): # Ex & Ez components - G.ID[0, :, 0, :] = 0 - G.ID[0, :, 1, :] = 0 - G.ID[2, :, 0, :] = 0 - G.ID[2, :, 1, :] = 0 + self.ID[0, :, 0, :] = 0 + self.ID[0, :, 1, :] = 0 + self.ID[2, :, 0, :] = 0 + self.ID[2, :, 1, :] = 0 - def tmz(): + def tmz(self): # Ex & Ey components - G.ID[0, :, :, 0] = 0 - G.ID[0, :, :, 1] = 0 - G.ID[1, :, :, 0] = 0 - G.ID[1, :, :, 1] = 0 + self.ID[0, :, :, 0] = 0 + self.ID[0, :, :, 1] = 0 + self.ID[1, :, :, 0] = 0 + self.ID[1, :, :, 1] = 0 def reset_fields(self): @@ -319,10 +316,8 @@ class FDTDGrid(Grid): self.dt = round_value(self.dt, decimalplaces=d.getcontext().prec - 1) - def dispersion_analysis(G): - """ - Analysis of numerical dispersion (Taflove et al, 2005, p112) - + """Analysis of numerical dispersion (Taflove et al, 2005, p112) - worse case of maximum frequency and minimum wavelength Args: @@ -449,12 +444,11 @@ def Ix(x, y, z, Hx, Hy, Hz, G): Args: x, y, z (float): Coordinates of position in grid. Hx, Hy, Hz (memory view): numpy array of magnetic field values. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ if y == 0 or z == 0: Ix = 0 - else: Ix = G.dy * (Hy[x, y, z - 1] - Hy[x, y, z]) + G.dz * (Hz[x, y, z] - Hz[x, y - 1, z]) @@ -467,12 +461,11 @@ def Iy(x, y, z, Hx, Hy, Hz, G): Args: x, y, z (float): Coordinates of position in grid. Hx, Hy, Hz (memory view): numpy array of magnetic field values. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ if x == 0 or z == 0: Iy = 0 - else: Iy = G.dx * (Hx[x, y, z] - Hx[x, y, z - 1]) + G.dz * (Hz[x - 1, y, z] - Hz[x, y, z]) @@ -485,12 +478,11 @@ def Iz(x, y, z, Hx, Hy, Hz, G): Args: x, y, z (float): Coordinates of position in grid. Hx, Hy, Hz (memory view): numpy array of magnetic field values. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ if x == 0 or y == 0: Iz = 0 - else: Iz = G.dx * (Hx[x, y - 1, z] - Hx[x, y, z]) + G.dy * (Hy[x, y, z] - Hy[x - 1, y, z]) diff --git a/gprMax/input_cmd_funcs.py b/gprMax/input_cmd_funcs.py deleted file mode 100644 index 6105bd1f..00000000 --- a/gprMax/input_cmd_funcs.py +++ /dev/null @@ -1,697 +0,0 @@ -# Copyright (C) 2015-2019: The University of Edinburgh -# Authors: Craig Warren and Antonis Giannopoulos -# -# This file is part of gprMax. -# -# gprMax is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# gprMax is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with gprMax. If not, see . - -import sys -from collections import namedtuple - -"""This module contains functional forms of some of the most commonly used gprMax -commands. It can be useful to use these within Python scripting in an input file. -For convenience, x, y, z coordinates are lumped in a namedtuple Coordinate: ->>> c = Coordinate(0.1, 0.2, 0.3) ->>> c -Coordinate(x=0.1, y=0.2, z=0.3) ->>> str(c) -'0.1 0.2 0.3' -# which can be accessed as a normal tuple: ->>> print c[0], c[1], c[2] -0.1 0.2 0.3 -# or more explicitly ->>> print c.x, c.y, c.z -0.1 0.2 0.3 -""" - -Coordinate_tuple = namedtuple('Coordinate', ['x', 'y', 'z']) - - -class Coordinate(Coordinate_tuple): - """Subclass of a namedtuple where __str__ outputs 'x y z'""" - - def __str__(self): - return '{:g} {:g} {:g}'.format(self.x, self.y, self.z) - - -def command(cmd, *parameters): - """ - Helper function. Prints the gprMax #: . None is ignored - in the output. - - Args: - cmd (str): the gprMax cmd string to be printed - *parameters: one or more strings as arguments, any None values are - ignored - - Returns: - s (str): the printed string - """ - - # remove Nones - filtered = filter(lambda x: x is not None, parameters) - # convert to str - filtered_str = map(str, filtered) - # convert to list - filtered_list = list(filtered_str) - try: - s = '#{}: {}'.format(cmd, " ".join(filtered_list)) - except TypeError as e: - # append info about cmd and parameters to the exception: - if not e.args: - e.args = ('', ) - additional_info = "Creating cmd = #{} with parameters {} -> {} failed".format(cmd, parameters, filtered_list) - e.args = e.args + (additional_info,) - raise e - # and now we can print it: - print(s) - return s - - -def rotate90_point(x, y, rotate90origin=()): - """Rotates a point 90 degrees CCW in the x-y plane. - - Args: - x, y (float): Coordinates. - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - xrot, yrot (float): Rotated coordinates. - """ - - # Translate point to origin - x -= rotate90origin[0] - y -= rotate90origin[1] - - # 90 degree CCW rotation and translate back - xrot = -y + rotate90origin[0] - yrot = x + rotate90origin[1] - - return xrot, yrot - - -def rotate90_edge(xs, ys, xf, yf, polarisation, rotate90origin): - """Rotates an edge or edge-like object/source 90 degrees CCW in the x-y plane. - - Args: - xs, ys, xf, yf (float): Start and finish coordinates. - polarisation (str): is the polarisation and can be 'x', 'y', or 'z'. - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - xs, ys, xf, yf (float): Rotated start and finish coordinates. - """ - - xsnew, ysnew = rotate90_point(xs, ys, rotate90origin) - xfnew, yfnew = rotate90_point(xf, yf, rotate90origin) - - # Swap coordinates for original y-directed edge, original x-directed - # edge does not require this. - if polarisation == 'y': - xs = xfnew - xf = xsnew - ys = ysnew - yf = yfnew - else: - xs = xsnew - xf = xfnew - ys = ysnew - yf = yfnew - - return xs, ys, xf, yf - - -def rotate90_plate(xs, ys, xf, yf, rotate90origin): - """Rotates an plate or plate-like object 90 degrees CCW in the x-y plane. - - Args: - xs, ys, xf, yf (float): Start and finish coordinates. - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - xs, ys, xf, yf (float): Rotated start and finish coordinates. - """ - - xsnew, ysnew = rotate90_point(xs, ys, rotate90origin) - xfnew, yfnew = rotate90_point(xf, yf, rotate90origin) - - # Swap x-coordinates to correctly specify plate - xs = xfnew - xf = xsnew - ys = ysnew - yf = yfnew - - return xs, ys, xf, yf - - -def domain(x, y, z): - """Prints the gprMax #domain command. - - Args: - x, y, z (float): Extent of the domain in the x, y, and z directions. - - Returns: - domain (Coordinate): Namedtuple of the extent of the domain. - """ - - domain = Coordinate(x, y, z) - command('domain', domain) - - return domain - - -def dx_dy_dz(x, y, z): - """Prints the gprMax #dx_dy_dz command. - - Args: - x, y, z (float): Spatial resolution in the x, y, and z directions. - - Returns: - dx_dy_dz (float): Tuple of the spatial resolutions. - """ - - dx_dy_dz = Coordinate(x, y, z) - command('dx_dy_dz', dx_dy_dz) - - return dx_dy_dz - - -def time_window(time_window): - """Prints the gprMax #time_window command. - - Args: - time_window (float): Duration of simulation. - - Returns: - time_window (float): Duration of simulation. - """ - - command('time_window', time_window) - - return time_window - - -def material(permittivity, conductivity, permeability, magconductivity, name): - """Prints the gprMax #material command. - - Args: - permittivity (float): Relative permittivity of the material. - conductivity (float): Conductivity of the material. - permeability (float): Relative permeability of the material. - magconductivity (float): Magnetic loss of the material. - name (str): Material identifier. - """ - - command('material', permittivity, conductivity, permeability, magconductivity, name) - - -def geometry_view(xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, type='n'): - """Prints the gprMax #geometry_view command. - - Args: - xs, ys, zs, xf, yf, zf (float): Start and finish coordinates. - dx, dy, dz (float): Spatial discretisation of geometry view. - filename (str): Filename where geometry file information will be stored. - type (str): 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). - - Returns: - s, f, d (tuple): 3 namedtuple Coordinate for the start, - finish coordinates and spatial discretisation - """ - - s = Coordinate(xs, ys, zs) - f = Coordinate(xf, yf, zf) - d = Coordinate(dx, dy, dz) - command('geometry_view', s, f, d, filename, type) - - return s, f, d - - -def snapshot(xs, ys, zs, xf, yf, zf, dx, dy, dz, time, filename): - """Prints the gprMax #snapshot command. - - Args: - xs, ys, zs, xf, yf, zf (float): Start and finish coordinates. - dx, dy, dz (float): Spatial discretisation of geometry view. - time (float): Time in seconds (float) or the iteration number - (integer) which denote the point in time at which the - snapshot will be taken. - filename (str): Filename where geometry file information will be stored. - - Returns: - s, f, d (tuple): 3 namedtuple Coordinate for the start, - finish coordinates and spatial discretisation - """ - - s = Coordinate(xs, ys, zs) - f = Coordinate(xf, yf, zf) - d = Coordinate(dx, dy, dz) - - if '.' in str(time) or 'e' in str(time): - time = '{:g}'.format(float(time)) - else: - time = '{:d}'.format(int(time)) - - command('snapshot', s, f, d, time, filename) - - return s, f, d - - -def edge(xs, ys, zs, xf, yf, zf, material, rotate90origin=()): - """Prints the gprMax #edge command. - - Args: - xs, ys, zs, xf, yf, zf (float): Start and finish coordinates. - material (str): Material identifier. - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - s, f (tuple): 2 namedtuple Coordinate for the start and finish coordinates - """ - - if rotate90origin: - if xs == xf: - polarisation = 'y' - else: - polarisation = 'x ' - xs, ys, xf, yf = rotate90_edge(xs, ys, xf, yf, polarisation, rotate90origin) - - s = Coordinate(xs, ys, zs) - f = Coordinate(xf, yf, zf) - command('edge', s, f, material) - - return s, f - - -def plate(xs, ys, zs, xf, yf, zf, material, rotate90origin=()): - """Prints the gprMax #plate command. - - Args: - xs, ys, zs, xf, yf, zf (float): Start and finish coordinates. - material (str): Material identifier(s). - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - s, f (tuple): 2 namedtuple Coordinate for the start and finish coordinates - """ - - if rotate90origin: - xs, ys, xf, yf = rotate90_plate(xs, ys, xf, yf, rotate90origin) - - s = Coordinate(xs, ys, zs) - f = Coordinate(xf, yf, zf) - command('plate', s, f, material) - - return s, f - - -def triangle(x1, y1, z1, x2, y2, z2, x3, y3, z3, thickness, material, averaging='', rotate90origin=()): - """Prints the gprMax #triangle command. - - Args: - x1, y1, z1, x2, y2, z2, x3, y3, z3 (float): Coordinates of the vertices. - thickness (float): Thickness for a triangular prism, or zero for a triangular patch. - material (str): Material identifier(s). - averaging (str): Turn averaging on or off. - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - v1, v2, v3 (tuple): 3 namedtuple Coordinate for the vertices - """ - - if rotate90origin: - x1, y1 = rotate90_point(x1, y1, rotate90origin) - x2, y2 = rotate90_point(x2, y2, rotate90origin) - x3, y3 = rotate90_point(x3, y3, rotate90origin) - - v1 = Coordinate(x1, y1, z1) - v2 = Coordinate(x2, y2, z2) - v3 = Coordinate(x3, y3, z3) - command('triangle', v1, v2, v3, thickness, material, averaging) - - return v1, v2, v3 - - -def box(xs, ys, zs, xf, yf, zf, material, averaging='', rotate90origin=()): - """Prints the gprMax #box command. - - Args: - xs, ys, zs, xf, yf, zf (float): Start and finish coordinates. - material (str): Material identifier(s). - averaging (str): Turn averaging on or off. - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - s, f (tuple): 2 namedtuple Coordinate for the start and finish coordinates - """ - - if rotate90origin: - xs, ys, xf, yf = rotate90_plate(xs, ys, xf, yf, rotate90origin) - - s = Coordinate(xs, ys, zs) - f = Coordinate(xf, yf, zf) - command('box', s, f, material, averaging) - - return s, f - - -def sphere(x, y, z, radius, material, averaging=''): - """Prints the gprMax #sphere command. - - Args: - x, y, z (float): Coordinates of the centre of the sphere. - radius (float): Radius. - material (str): Material identifier(s). - averaging (str): Turn averaging on or off. - - Returns: - c (tuple): namedtuple Coordinate for the center of the sphere - """ - - c = Coordinate(x, y, z) - command('sphere', c, radius, material, averaging) - - return c - - -def cylinder(x1, y1, z1, x2, y2, z2, radius, material, averaging='', rotate90origin=()): - """Prints the gprMax #cylinder command. - - Args: - x1, y1, z1, x2, y2, z2 (float): Coordinates of the centres of the two faces of the cylinder. - radius (float): Radius. - material (str): Material identifier(s). - averaging (str): Turn averaging on or off. - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - c1, c2 (tuple): 2 namedtuple Coordinate for the centres of the two faces of the cylinder. - """ - - if rotate90origin: - x1, y1 = rotate90_point(x1, y1, rotate90origin) - x2, y2 = rotate90_point(x2, y2, rotate90origin) - - c1 = Coordinate(x1, y1, z1) - c2 = Coordinate(x2, y2, z2) - command('cylinder', c1, c2, radius, material, averaging) - - return c1, c2 - - -def cylindrical_sector(axis, ctr1, ctr2, t1, t2, radius, - startingangle, sweptangle, material, averaging=''): - """Prints the gprMax #cylindrical_sector command. - - Args: - axis (str): Axis of the cylinder from which the sector is defined and - can be x, y, or z. - ctr1, ctr2 (float): Coordinates of the centre of the cylindrical sector. - t1, t2 (float): 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). - radius (float): Radius. - startingangle (float): Starting angle (in degrees) for the cylindrical - sector (with zero degrees defined on the positive first axis of - the plane of the cylindrical sector). - sweptangle (float): Angle (in degrees) swept by the cylindrical sector - (the finishing angle of the sector is always anti-clockwise - from the starting angle). - material (str): Material identifier(s). - averaging (str): Turn averaging on or off. - """ - - command('cylindrical_sector', axis, ctr1, ctr2, t1, t2, radius, startingangle, sweptangle, material, averaging) - - -def excitation_file(file1): - """Prints the #excitation_file: command. - - Args: - file1 (str): filename - - Returns: - file1 (str): filename - """ - - command('excitation_file', file1) - - return file1 - - -def waveform(shape, amplitude, frequency, identifier): - """Prints the #waveform: shape amplitude frequency identifier - - Args: - shape (str): is the type of waveform - amplitude (float): is the amplitude of the waveform. - frequency (float): is the frequency of the waveform in Hertz. - identifier (str): is an identifier for the waveform used to assign it to a source. - - Returns: - identifier (str): is an identifier for the waveform used to assign it to a source. - """ - - command('waveform', shape, amplitude, frequency, identifier) - - return identifier - - -def hertzian_dipole(polarisation, f1, f2, f3, identifier, - t0=None, t_remove=None, dxdy=None, rotate90origin=()): - """Prints the #hertzian_dipole: polarisation, f1, f2, f3, identifier, [t0, t_remove] - - Args: - polarisation (str): is the polarisation of the source and can be 'x', 'y', or 'z'. - f1 f2 f3 (float): are the coordinates (x,y,z) of the source in the model. - identifier (str): is the identifier of the waveform that should be used with the source. - t0 (float): is an optinal argument for the time delay in starting the source. - t_remove (float): is a time to remove the source. - dxdy (float): Tuple of x-y spatial resolutions. For rotation purposes only. - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - coordinates (tuple): namedtuple Coordinate of the source location - """ - - if rotate90origin: - if polarisation == 'x': - xf = f1 + dxdy[0] - yf = f2 - newpolarisation = 'y' - elif polarisation == 'y': - xf = f1 - yf = f2 + dxdy[1] - newpolarisation = 'x' - - f1, f2, xf, yf = rotate90_edge(f1, f2, xf, yf, polarisation, rotate90origin) - polarisation = newpolarisation - - c = Coordinate(f1, f2, f3) - # since command ignores None, this is safe: - command('hertzian_dipole', polarisation, str(c), identifier, t0, t_remove) - - return c - - -def magnetic_dipole(polarisation, f1, f2, f3, identifier, - t0=None, t_remove=None, dxdy=None, rotate90origin=()): - """Prints the #magnetic_dipole: polarisation, f1, f2, f3, identifier, [t0, t_remove] - - Args: - polarisation (str): is the polarisation of the source and can be 'x', 'y', or 'z'. - f1 f2 f3 (float): are the coordinates (x,y,z) of the source in the model. - identifier (str): is the identifier of the waveform that should be used with the source. - t0 (float): is an optinal argument for the time delay in starting the source. - t_remove (float): is a time to remove the source. - dxdy (float): Tuple of x-y spatial resolutions. For rotation purposes only. - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - coordinates (tuple): namedtuple Coordinate of the source location - """ - - if rotate90origin: - if polarisation == 'x': - xf = f1 + dxdy[0] - yf = f2 - newpolarisation = 'y' - elif polarisation == 'y': - xf = f1 - yf = f2 + dxdy[1] - newpolarisation = 'x' - - f1, f2, xf, yf = rotate90_edge(f1, f2, xf, yf, polarisation, rotate90origin) - polarisation = newpolarisation - - c = Coordinate(f1, f2, f3) - # since command ignores None, this is safe: - command('magnetic_dipole', polarisation, str(c), identifier, t0, t_remove) - - return c - - -def voltage_source(polarisation, f1, f2, f3, resistance, identifier, - t0=None, t_remove=None, dxdy=None, rotate90origin=()): - """Prints the #voltage_source: polarisation, f1, f2, f3, resistance, identifier, [t0, t_remove] - - Args: - polarisation (str): is the polarisation of the source and can be 'x', 'y', or 'z'. - f1 f2 f3 (float): are the coordinates (x,y,z) of the source in the model. - identifier (str): is the identifier of the waveform that should be used with the source. - resistance (float): is the internal resistance of the voltage source. - t0 (float): is an optinal argument for the time delay in starting the source. - t_remove (float): is a time to remove the source. - dxdy (float): Tuple of x-y spatial resolutions. For rotation purposes only. - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - coordinates (tuple): namedtuple Coordinate of the source location - """ - - if rotate90origin: - if polarisation == 'x': - xf = f1 + dxdy[0] - yf = f2 - newpolarisation = 'y' - elif polarisation == 'y': - xf = f1 - yf = f2 + dxdy[1] - newpolarisation = 'x' - - f1, f2, xf, yf = rotate90_edge(f1, f2, xf, yf, polarisation, rotate90origin) - polarisation = newpolarisation - - c = Coordinate(f1, f2, f3) - # since command ignores None, this is safe: - command('voltage_source', polarisation, str(c), resistance, identifier, t0, t_remove) - - return c - - -def transmission_line(polarisation, f1, f2, f3, resistance, identifier, - t0=None, t_remove=None, dxdy=None, rotate90origin=()): - """Prints the #transmission_line: polarisation, f1, f2, f3, resistance, identifier, [t0, t_remove] - - Args: - polarisation (str): is the polarisation of the source and can be 'x', 'y', or 'z'. - f1 f2 f3 (float): are the coordinates (x,y,z) of the source in the model. - identifier (str): is the identifier of the waveform that should be used with the source. - resistance (float): is the characteristic resistance of the transmission_line. - t0 (float): is an optinal argument for the time delay in starting the source. - t_remove (float): is a time to remove the source. - dxdy (float): Tuple of x-y spatial resolutions. For rotation purposes only. - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - coordinates (tuple): namedtuple Coordinate of the source location - """ - - if rotate90origin: - if polarisation == 'x': - xf = f1 + dxdy[0] - yf = f2 - newpolarisation = 'y' - elif polarisation == 'y': - xf = f1 - yf = f2 + dxdy[1] - newpolarisation = 'x' - - f1, f2, xf, yf = rotate90_edge(f1, f2, xf, yf, polarisation, rotate90origin) - polarisation = newpolarisation - - c = Coordinate(f1, f2, f3) - # since command ignores None, this is safe: - command('transmission_line', polarisation, str(c), resistance, identifier, t0, t_remove) - - return c - - -def rx(x, y, z, identifier=None, to_save=None, polarisation=None, dxdy=None, rotate90origin=()): - """Prints the #rx: x, y, z, [identifier, to_save] command. - - Args: - x, y, z (float): are the coordinates (x,y,z) of the receiver in the model. - identifier (str): is the optional identifier of the receiver - to_save (list): is a list of outputs with this receiver. It can be - any selection from 'Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', 'Ix', 'Iy', or 'Iz'. - polarisation (str): is the polarisation of the source - and can be 'x', 'y', or 'z'. For rotation purposes only. - dxdy (float): Tuple of x-y spatial resolutions. For rotation purposes only. - rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane. - - Returns: - coordinates (tuple): namedtuple Coordinate of the receiver location - """ - - if rotate90origin: - if polarisation == 'x': - try: - xf = x + dxdy[0] - except Exception as e: - raise ValueError('With polarization = x, a dxdy[0] float \ - values is required, got dxdy=%s' % dxdy) from e - yf = y - elif polarisation == 'y': - xf = x - try: - yf = y + dxdy[1] - except Exception as e: - raise ValueError('With polarization = y, a dxdy[1] float \ - values is required, got dxdy=%s' % dxdy) from e - - x, y, xf, yf = rotate90_edge(x, y, xf, yf, polarisation, rotate90origin) - - c = Coordinate(x, y, z) - to_save_str = '' - if to_save is not None: - to_save_str = ''.join(to_save) - - command('rx', str(c), identifier, to_save_str) - - return c - - -def src_steps(dx=0, dy=0, dz=0): - """Prints the #src_steps: dx, dy, dz command. - - Args: - dx, dy, dz (float): are the increments in (x, y, z) to - move all simple sources or all receivers. - - Returns: - coordinates (tuple): namedtuple Coordinate of the increments - """ - - c = Coordinate(dx, dy, dz) - command('src_steps', str(c)) - - return c - - -def rx_steps(dx=0, dy=0, dz=0): - """Prints the #rx_steps: dx, dy, dz command. - - Args: - dx, dy, dz (float): are the increments in (x, y, z) to move all simple sources or all receivers. - - Returns: - coordinates (tuple): namedtuple Coordinate of the increments - """ - - c = Coordinate(dx, dy, dz) - command('rx_steps', str(c)) - return c diff --git a/gprMax/input_cmds_file.py b/gprMax/input_cmds_file.py index 32cfec75..e2e32ed4 100644 --- a/gprMax/input_cmds_file.py +++ b/gprMax/input_cmds_file.py @@ -22,7 +22,6 @@ from io import StringIO import gprMax.config as config from .exceptions import CmdInputError - from .input_cmds_geometry import process_geometrycmds from .input_cmds_multiuse import process_multicmds from .input_cmds_singleuse import process_singlecmds @@ -113,8 +112,7 @@ def process_python_include_code(inputfile, usernamespace): def process_include_files(hashcmds, inputfile): - """ - Looks for and processes any include file commands and insert + """Looks for and processes any include file commands and insert the contents of the included file at that location. Args: @@ -157,9 +155,8 @@ def process_include_files(hashcmds, inputfile): def write_processed_file(processedlines, appendmodelnumber): - """ - Writes an input file after any Python code and include commands - in the original input file have been processed. + """Writes an input file after any Python code and include commands + in the original input file have been processed. Args: processedlines (list): Input commands after after processing any @@ -177,8 +174,7 @@ def write_processed_file(processedlines, appendmodelnumber): def check_cmd_names(processedlines, checkessential=True): - """ - Checks the validity of commands, i.e. are they gprMax commands, + """Checks the validity of commands, i.e. are they gprMax commands, and that all essential commands are present. Args: @@ -253,9 +249,17 @@ def check_cmd_names(processedlines, checkessential=True): def get_user_objects(processedlines, check=True): + """Make a list of all user objects. - # Check validity of command names and that - # essential commands are present + Args: + processedlines (list): Input commands after Python processing. + check (bool): Whether to check for essential commands or not. + + Returns: + user_objs (list): All user objects. + """ + + # Check validity of command names and that essential commands are present parsed_commands = check_cmd_names(processedlines, checkessential=check) # Process parameters for commands that can only occur once in the model @@ -274,9 +278,16 @@ def get_user_objects(processedlines, check=True): def parse_hash_commands(model_config, G, scene): - """Parse user hash commands and add them to the scene.""" - # Add the current model run to namespace that can be accessed by - # user in any Python code blocks in input file + """Parse user hash commands and add them to the scene. + + Args: + model_config (ModelConfig): Model level configuration object. + G (Grid): Holds essential parameters describing the model. + scene (Scene): Scene object. + + Returns: + scene (Scene): Scene object. + """ sim_config = model_config.sim_config @@ -292,8 +303,7 @@ def parse_hash_commands(model_config, G, scene): for key, value in sorted(usernamespace.items()): if key != '__builtins__': uservars += '{}: {}, '.format(key, value) - usv_s = """Constants/variables used/available for Python \ - scripting: {{{}}}\n""" + usv_s = """Constants/variables used/available for Python scripting: {{{}}}\n""" print(usv_s.format(uservars[:-2])) # Write a file containing the input commands after Python or include @@ -309,8 +319,8 @@ def parse_hash_commands(model_config, G, scene): class Capturing(list): + """Context manager to capture standard output stream.""" # https://stackoverflow.com/questions/16571150/how-to-capture-stdout-output-from-a-python-function-call - """Context manager to capture standard output stream""" def __enter__(self): self._stdout = sys.stdout @@ -325,10 +335,12 @@ class Capturing(list): def user_libs_fn_to_scene_obj(f, *args, **kwargs): """Function to convert library functions in the user_libs directory - into geometry objects which can be added to the scene""" + into geometry objects which can be added to the scene. + """ with Capturing() as str_cmds: f(*args, **kwargs) user_objects = get_user_objects(str_cmds, check=False) + return user_objects diff --git a/gprMax/input_cmds_geometry.py b/gprMax/input_cmds_geometry.py index a3793082..74e8ca36 100644 --- a/gprMax/input_cmds_geometry.py +++ b/gprMax/input_cmds_geometry.py @@ -15,9 +15,9 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -import sys -from .utilities import get_terminal_width +import sys +import numpy as np from .cmds_geometry.geometry_objects_read import GeometryObjectsRead from .cmds_geometry.edge import Edge @@ -31,14 +31,9 @@ from .cmds_geometry.sphere import Sphere from .cmds_geometry.add_surface_roughness import AddSurfaceRoughness from .cmds_geometry.add_surface_water import AddSurfaceWater from .cmds_geometry.add_grass import AddGrass - from .utilities import round_value -from tqdm import tqdm -import numpy as np - - def process_geometrycmds(geometry): """ This function checks the validity of command parameters, creates instances @@ -46,21 +41,15 @@ def process_geometrycmds(geometry): solid, rigid and ID. Args: - geometry (list): Geometry commands in the model + geometry (list): Geometry commands in the model, + + Returns: + scene_objects (list): Holds objects in scene. """ scene_objects = [] - # Disable progress bar if on Windows as it does not update properly - # when messages are printed - #if sys.platform == 'win32': -# tqdmdisable = True - #else: - # tqdmdisable = G.tqdmdisable - - tqdmdisable = False - - for object in tqdm(geometry, desc='Processing geometry related cmds', unit='cmds', ncols=get_terminal_width() - 1, file=sys.stdout, disable=tqdmdisable): + for object in geometry: tmp = object.split() if tmp[0] == '#geometry_objects_read:': diff --git a/gprMax/input_cmds_multiuse.py b/gprMax/input_cmds_multiuse.py index 164b3549..7098e672 100644 --- a/gprMax/input_cmds_multiuse.py +++ b/gprMax/input_cmds_multiuse.py @@ -32,7 +32,7 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -from .exceptions import CmdInputError + from .cmds_multiple import Waveform from .cmds_multiple import VoltageSource from .cmds_multiple import HertzianDipole @@ -48,6 +48,7 @@ from .cmds_multiple import GeometryView from .cmds_multiple import GeometryObjectsWrite from .cmds_multiple import PMLCFS from .cmds_multiple import Rx +from .exceptions import CmdInputError def process_multicmds(multicmds): @@ -57,7 +58,9 @@ def process_multicmds(multicmds): Args: multicmds (dict): Commands that can have multiple instances in the model. - G (class): Grid class instance - holds essential parameters describing the model. + + Returns: + scene_objects (list): Holds objects in scene. """ scene_objects = [] @@ -333,6 +336,7 @@ def process_multicmds(multicmds): return scene_objects + def process_subgrid_hsg(cmdinstance): pass diff --git a/gprMax/input_cmds_singleuse.py b/gprMax/input_cmds_singleuse.py index 6fdf3787..a1efa76a 100644 --- a/gprMax/input_cmds_singleuse.py +++ b/gprMax/input_cmds_singleuse.py @@ -32,6 +32,7 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . + from .cmds_single_use import Messages from .cmds_single_use import Title from .cmds_single_use import NumThreads @@ -48,16 +49,21 @@ from .exceptions import CmdInputError def process_singlecmds(singlecmds): - """Checks the validity of command parameters and creates instances of classes of parameters. + """ + Checks the validity of command parameters and creates instances of + classes of parameters. Args: singlecmds (dict): Commands that can only occur once in the model. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. + + Returns: + scene_objects (list): Holds objects in scene. """ + scene_objects = [] # Check validity of command parameters in order needed - # messages cmd = '#messages' if singlecmds[cmd] is not None: tmp = singlecmds[cmd].split() @@ -67,12 +73,18 @@ def process_singlecmds(singlecmds): messages = Messages(yn=str(tmp[0])) scene_objects.append(messages) - # Title cmd = '#title' if singlecmds[cmd] is not None: title = Title(name=str(singlecmds[cmd])) scene_objects.append(title) + # Set the output directory + cmd = '#output_dir' + if singlecmds[cmd] is not None: + output_dir = OutputDir(dir=singlecmds[cmd]) + scene_objects.append(output_dir) + + # Number of threads for CPU-based (OpenMP) parallelised parts of code cmd = '#num_threads' if singlecmds[cmd] is not None: tmp = tuple(int(x) for x in singlecmds[cmd].split()) @@ -93,7 +105,6 @@ def process_singlecmds(singlecmds): discretisation = Discretisation(p1=dl) scene_objects.append(discretisation) - # Domain cmd = '#domain' if singlecmds[cmd] is not None: tmp = [float(x) for x in singlecmds[cmd].split()] @@ -104,14 +115,12 @@ def process_singlecmds(singlecmds): domain = Domain(p1=p1) scene_objects.append(domain) - # Time step stability factor cmd = '#time_step_stability_factor' if singlecmds[cmd] is not None: tmp = tuple(float(x) for x in singlecmds[cmd].split()) tsf = TimeStepStabilityFactor(f=tmp[0]) scene_objects.append(tsf) - # Time window cmd = '#time_window' if singlecmds[cmd] is not None: tmp = singlecmds[cmd].split() @@ -132,25 +141,23 @@ def process_singlecmds(singlecmds): scene_objects.append(tw) - # PML - cmd = '#pml_cells' - if singlecmds[cmd] is not None: - tmp = singlecmds[cmd].split() - if len(tmp) != 1 and len(tmp) != 6: - raise CmdInputError(cmd + ' requires either one or six parameter(s)') - if len(tmp) == 1: - pml_cells = PMLCells(thickness=int(tmp[0])) - else: - pml_cells = PMLCells(x0=int(tmp[0]), - y0=int(tmp[1]), - z0=int(tmp[2]), - xmax=int(tmp[3]), - ymax=int(tmp[4]), - zmax=int(tmp[5])) + cmd = '#pml_cells' + if singlecmds[cmd] is not None: + tmp = singlecmds[cmd].split() + if len(tmp) != 1 and len(tmp) != 6: + raise CmdInputError(cmd + ' requires either one or six parameter(s)') + if len(tmp) == 1: + pml_cells = PMLCells(thickness=int(tmp[0])) + else: + pml_cells = PMLCells(x0=int(tmp[0]), + y0=int(tmp[1]), + z0=int(tmp[2]), + xmax=int(tmp[3]), + ymax=int(tmp[4]), + zmax=int(tmp[5])) - scene_objects.append(pml_cells) + scene_objects.append(pml_cells) - # src_steps cmd = '#src_steps' if singlecmds[cmd] is not None: tmp = singlecmds[cmd].split() @@ -160,7 +167,6 @@ def process_singlecmds(singlecmds): src_steps = SrcSteps(p1=p1) scene_objects.append(src_steps) - # rx_steps cmd = '#rx_steps' if singlecmds[cmd] is not None: tmp = singlecmds[cmd].split() @@ -185,10 +191,4 @@ def process_singlecmds(singlecmds): scene_objects.append(ex_file) - # Set the output directory - cmd = '#output_dir' - if singlecmds[cmd] is not None: - output_dir = OutputDir(dir=singlecmds[cmd]) - scene_objects.append(output_dir) - return scene_objects diff --git a/gprMax/materials.py b/gprMax/materials.py index 15558e5e..d75f22b1 100644 --- a/gprMax/materials.py +++ b/gprMax/materials.py @@ -49,8 +49,7 @@ class Material(object): """Calculates the magnetic update coefficients of the material. Args: - G (class): Grid class instance - holds essential parameters - describing the model. + G (Grid): Holds essential parameters describing the model. """ HA = (config.m0 * self.mr / G.dt) + 0.5 * self.sm @@ -65,8 +64,7 @@ class Material(object): """Calculates the electric update coefficients of the material. Args: - G (class): Grid class instance - holds essential parameters - describing the model. + G (Grid): Holds essential parameters describing the model. """ EA = (config.e0 * self.er / G.dt) + 0.5 * self.se @@ -128,8 +126,7 @@ class DispersiveMaterial(Material): """Calculates the electric update coefficients of the material. Args: - G (class): Grid class instance - holds essential parameters - describing the model. + G (Grid): Holds essential parameters describing the model. """ # The implementation of the dispersive material modelling comes from the @@ -211,8 +208,7 @@ def process_materials(G): store in arrays, and build text list of materials/properties Args: - G (class): Grid class instance - holds essential parameters - describing the model. + G (Grid): Holds essential parameters describing the model. Returns: materialsdata (list): List of material IDs, names, and properties to @@ -311,7 +307,7 @@ class PeplinskiSoil(object): Args: nbins (int): Number of bins to use to create the different materials. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ # Debye model properties of water @@ -373,15 +369,22 @@ class PeplinskiSoil(object): muiter.iternext() + def create_built_in_materials(G): + """Create pre-defined (built-in) materials. + + Args: + G (Grid): Holds essential parameters describing the model. + """ + G.n_built_in_materials = len(G.materials) - # Create built-in materials m = Material(0, 'pec') m.se = float('inf') m.type = 'builtin' m.averagable = False G.materials.append(m) + m = Material(1, 'free_space') m.type = 'builtin' G.materials.append(m) diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py index 0fb996f9..8cc65f23 100644 --- a/gprMax/model_build_run.py +++ b/gprMax/model_build_run.py @@ -60,9 +60,11 @@ from .pml import pml_information from .receivers import gpu_initialise_rx_arrays from .receivers import gpu_get_rx_array from .receivers import Rx +from .scene import Scene from .snapshots import Snapshot from .snapshots import gpu_initialise_snapshot_array from .snapshots import gpu_get_snapshot_array +from .solvers import create_solver from .sources import gpu_initialise_src_arrays from .utilities import get_terminal_width from .utilities import human_size @@ -70,12 +72,10 @@ from .utilities import open_path_file from .utilities import round32 from .utilities import timer from .utilities import Printer -from .scene import Scene -from .solvers import create_solver - class ModelBuildRun: + """Builds and runs (solves) a model.""" def __init__(self, G, sim_config, model_config): self.G = G @@ -86,20 +86,8 @@ class ModelBuildRun: self.p = None def build(self): - """Runs a model - processes the input file; builds the Yee cells; calculates update coefficients; runs main FDTD loop. + """Builds the Yee cells for a model.""" - Args: - args (dict): Namespace with command line arguments - currentmodelrun (int): Current model run number. - modelend (int): Number of last model to run. - numbermodelruns (int): Total number of model runs. - inputfile (object): File object for the input file. - usernamespace (dict): Namespace that can be accessed by user - in any Python code blocks in input file. - - Returns: - tsolve (int): Length of time (seconds) of main FDTD calculations - """ # Monitor memory usage self.p = psutil.Process() @@ -114,20 +102,20 @@ class ModelBuildRun: # Adjust position of simple sources and receivers if required if G.srcsteps[0] != 0 or G.srcsteps[1] != 0 or G.srcsteps[2] != 0: for source in itertools.chain(G.hertziandipoles, G.magneticdipoles): - if currentmodelrun == 1: - if source.xcoord + G.srcsteps[0] * modelend < 0 or source.xcoord + G.srcsteps[0] * modelend > G.nx or source.ycoord + G.srcsteps[1] * modelend < 0 or source.ycoord + G.srcsteps[1] * modelend > G.ny or source.zcoord + G.srcsteps[2] * modelend < 0 or source.zcoord + G.srcsteps[2] * modelend > G.nz: + if self.model_config.i == 0: + if source.xcoord + G.srcsteps[0] * self.sim_config.model_end < 0 or source.xcoord + G.srcsteps[0] * self.sim_config.model_end > G.nx or source.ycoord + G.srcsteps[1] * self.sim_config.model_end < 0 or source.ycoord + G.srcsteps[1] * self.sim_config.model_end > G.ny or source.zcoord + G.srcsteps[2] * self.sim_config.model_end < 0 or source.zcoord + G.srcsteps[2] * self.sim_config.model_end > G.nz: raise GeneralError('Source(s) will be stepped to a position outside the domain.') - source.xcoord = source.xcoordorigin + (currentmodelrun - 1) * G.srcsteps[0] - source.ycoord = source.ycoordorigin + (currentmodelrun - 1) * G.srcsteps[1] - source.zcoord = source.zcoordorigin + (currentmodelrun - 1) * G.srcsteps[2] + source.xcoord = source.xcoordorigin + (self.model_config.i) * G.srcsteps[0] + source.ycoord = source.ycoordorigin + (self.model_config.i) * G.srcsteps[1] + source.zcoord = source.zcoordorigin + (self.model_config.i) * G.srcsteps[2] if G.rxsteps[0] != 0 or G.rxsteps[1] != 0 or G.rxsteps[2] != 0: for receiver in G.rxs: - if currentmodelrun == 1: - if receiver.xcoord + G.rxsteps[0] * modelend < 0 or receiver.xcoord + G.rxsteps[0] * modelend > G.nx or receiver.ycoord + G.rxsteps[1] * modelend < 0 or receiver.ycoord + G.rxsteps[1] * modelend > G.ny or receiver.zcoord + G.rxsteps[2] * modelend < 0 or receiver.zcoord + G.rxsteps[2] * modelend > G.nz: + if self.model_config.i == 0: + if receiver.xcoord + G.rxsteps[0] * self.sim_config.model_end < 0 or receiver.xcoord + G.rxsteps[0] * self.sim_config.model_end > G.nx or receiver.ycoord + G.rxsteps[1] * self.sim_config.model_end < 0 or receiver.ycoord + G.rxsteps[1] * self.sim_config.model_end > G.ny or receiver.zcoord + G.rxsteps[2] * self.sim_config.model_end < 0 or receiver.zcoord + G.rxsteps[2] * self.sim_config.model_end > G.nz: raise GeneralError('Receiver(s) will be stepped to a position outside the domain.') - receiver.xcoord = receiver.xcoordorigin + (currentmodelrun - 1) * G.rxsteps[0] - receiver.ycoord = receiver.ycoordorigin + (currentmodelrun - 1) * G.rxsteps[1] - receiver.zcoord = receiver.zcoordorigin + (currentmodelrun - 1) * G.rxsteps[2] + receiver.xcoord = receiver.xcoordorigin + (self.model_config.i) * G.rxsteps[0] + receiver.ycoord = receiver.ycoordorigin + (self.model_config.i) * G.rxsteps[1] + receiver.zcoord = receiver.zcoordorigin + (self.model_config.i) * G.rxsteps[2] # Write files for any geometry views and geometry object outputs if not (G.geometryviews or G.geometryobjectswrite) and self.sim_config.geometry_only and config.is_messages(): @@ -142,15 +130,15 @@ class ModelBuildRun: pbar = tqdm(total=geometryobject.datawritesize, unit='byte', unit_scale=True, desc='Writing geometry object file {}/{}, {}'.format(i + 1, len(G.geometryobjectswrite), os.path.split(geometryobject.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) geometryobject.write_hdf5(G, pbar) pbar.close() + if config.is_messages(): print() - # If only writing geometry information - if self.sim_config.geometry_only: - tsolve = 0 + # # If only writing geometry information + # if self.sim_config.geometry_only: + # tsolve = 0 def build_geometry(self): model_config = self.model_config sim_config = self.sim_config - G = self.G printer = Printer(config) @@ -158,7 +146,7 @@ class ModelBuildRun: scene = self.build_scene() - # combine available grids + # Combine available grids grids = [G] + G.subgrids gridbuilders = [GridBuilder(grid, self.printer) for grid in grids] @@ -199,6 +187,7 @@ class ModelBuildRun: printer.print('\nMemory (RAM) required - updated (snapshots): ~{}\n'.format(human_size(G.memoryusage))) + # Build materials for gb in gridbuilders: gb.build_materials() @@ -213,31 +202,32 @@ class ModelBuildRun: elif results['deltavp']: printer.print("\nNumerical dispersion analysis: estimated largest physical phase-velocity error is {:.2f}% in material '{}' whose wavelength sampled by {} cells. Maximum significant frequency estimated as {:g}Hz".format(results['deltavp'], results['material'].ID, results['N'], results['maxfreq'])) - def reuse_geometry(self): G = self.G - inputfilestr = '\n--- Model {}/{}, input file (not re-processed, i.e. geometry fixed): {}'.format(self.model_config.appendmodelnumber, self.sim_config.model_end, self.sim_config.input_file_path) - self.printer.print(Fore.GREEN + '{} {}\n'.format(self.model_config.inputfilestr, '-' * (get_terminal_width() - 1 - len(inputfilestr))) + Style.RESET_ALL) + # Reset iteration number + G.iteration = 0 + self.model_config.inputfilestr = '\n--- Model {}/{}, input file (not re-processed, i.e. geometry fixed): {}'.format(self.model_config.appendmodelnumber, self.sim_config.model_end, self.sim_config.input_file_path) + self.printer.print(Fore.GREEN + '{} {}'.format(self.model_config.inputfilestr, '-' * (get_terminal_width() - 1 - len(self.model_config.inputfilestr))) + Style.RESET_ALL) for grid in [G] + G.subgrids: grid.reset_fields() def build_scene(self): G = self.G - # api for multiple scenes / model runs + # API for multiple scenes / model runs scene = self.model_config.get_scene() - # if there is no scene - process the hash commands instead + # If there is no scene - process the hash commands instead if not scene: scene = Scene() - # parse the input file into user objects and add them to the scene + # Parse the input file into user objects and add them to the scene scene = parse_hash_commands(self.model_config, G, scene) # Creates the internal simulation objects. scene.create_internal_objects(G) + return scene def create_output_directory(self): - if self.G.outputdirectory: # Check and set output directory and filename try: @@ -245,17 +235,15 @@ class ModelBuildRun: self.printer.print('\nCreated output directory: {}'.format(self.G.outputdirectory)) except FileExistsError: pass - # modify the output path (hack) + # Modify the output path (hack) self.model_config.output_file_path_ext = Path(self.G.outputdirectory, self.model_config.output_file_path_ext) - def run_model(self, solver): - G = self.G - self.create_output_directory() - self.printer.print('\nOutput file: {}\n'.format(self.model_config.output_file_path_ext)) + self.printer.print('Output file: {}\n'.format(self.model_config.output_file_path_ext)) + # Run solver tsolve = self.solve(solver) # Write an output file in HDF5 format @@ -292,13 +280,12 @@ class ModelBuildRun: electric and magnetic field updates, and PML updates. Args: - currentmodelrun (int): Current model run number. - modelend (int): Number of last model to run. - G (class): Grid class instance - holds essential parameters describing the model. + solver (Solver): solver object. Returns: - tsolve (float): Time taken to execute solving (seconds) + tsolve (float): time taken to execute solving (seconds). """ + G = self.G if config.is_messages(): @@ -318,21 +305,21 @@ class GridBuilder: self.printer = printer def build_pmls(self): - grid = self.grid - # build the PMLS - pbar = tqdm(total=sum(1 for value in grid.pmlthickness.values() if value > 0), desc='Building {} Grid PML boundaries'.format(self.grid.name), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) + # Build the PMLS + pbar = tqdm(total=sum(1 for value in grid.pmlthickness.values() if value > 0), desc='Building {} Grid PML boundaries'.format(self.grid.name), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) for pml_id, thickness in grid.pmlthickness.items(): - build_pml(grid, pml_id, thickness) - pbar.update() + if thickness > 0: + build_pml(grid, pml_id, thickness) + pbar.update() pbar.close() def build_components(self): # Build the model, i.e. set the material properties (ID) for every edge # of every Yee cell self.printer.print('') - pbar = tqdm(total=2, desc='Building {} grid'.format(self.grid.name), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) + pbar = tqdm(total=2, desc='Building {} Grid'.format(self.grid.name), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars']) build_electric_components(self.grid.solid, self.grid.rigidE, self.grid.ID, self.grid) pbar.update() build_magnetic_components(self.grid.solid, self.grid.rigidH, self.grid.ID, self.grid) diff --git a/gprMax/pml.py b/gprMax/pml.py index 59d1a942..4484db73 100644 --- a/gprMax/pml.py +++ b/gprMax/pml.py @@ -177,8 +177,7 @@ class PML(object): def __init__(self, G, ID=None, direction=None, xs=0, xf=0, ys=0, yf=0, zs=0, zf=0): """ Args: - G (class): Grid class instance - holds essential parameters - describing the model. + G (Grid): Holds essential parameters describing the model. ID (str): Identifier for PML slab. direction (str): Direction of increasing absorption. xs, xf, ys, yf, zs, zf (float): Extent of the PML slab. @@ -249,8 +248,7 @@ class PML(object): Args: er (float): Average permittivity of underlying material mr (float): Average permeability of underlying material - G (class): Grid class instance - holds essential parameters - describing the model. + G (Grid): Holds essential parameters describing the model. """ self.ERA = np.zeros((len(self.CFS), self.thickness), @@ -313,7 +311,7 @@ class PML(object): """This functions updates electric field components with the PML correction. Args: - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ pmlmodule = 'gprMax.cython.pml_updates_electric_' + G.pmlformulation @@ -327,7 +325,7 @@ class PML(object): """This functions updates magnetic field components with the PML correction. Args: - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ pmlmodule = 'gprMax.cython.pml_updates_magnetic_' + G.pmlformulation @@ -341,7 +339,7 @@ class PML(object): """Set the blocks per grid size used for updating the PML field arrays on a GPU. Args: - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ self.bpg = (int(np.ceil(((self.EPhi1_gpu.shape[1] + 1) * (self.EPhi1_gpu.shape[2] + 1) * (self.EPhi1_gpu.shape[3] + 1)) / G.tpb[0])), 1, 1) @@ -394,8 +392,7 @@ class PML(object): correction on the GPU. Args: - G (class): Grid class instance - holds essential parameters - describing the model. + G (Grid): Holds essential parameters describing the model. """ self.update_electric_gpu(np.int32(self.xs), np.int32(self.xf), np.int32(self.ys), np.int32(self.yf), np.int32(self.zs), np.int32(self.zf), np.int32(self.EPhi1_gpu.shape[1]), np.int32(self.EPhi1_gpu.shape[2]), np.int32(self.EPhi1_gpu.shape[3]), np.int32(self.EPhi2_gpu.shape[1]), np.int32(self.EPhi2_gpu.shape[2]), np.int32(self.EPhi2_gpu.shape[3]), np.int32(self.thickness), G.ID_gpu.gpudata, G.Ex_gpu.gpudata, G.Ey_gpu.gpudata, G.Ez_gpu.gpudata, G.Hx_gpu.gpudata, G.Hy_gpu.gpudata, G.Hz_gpu.gpudata, self.EPhi1_gpu.gpudata, self.EPhi2_gpu.gpudata, self.ERA_gpu.gpudata, self.ERB_gpu.gpudata, self.ERE_gpu.gpudata, self.ERF_gpu.gpudata, floattype(self.d), block=G.tpb, grid=self.bpg) @@ -405,14 +402,13 @@ class PML(object): correction on the GPU. Args: - G (class): Grid class instance - holds essential parameters - describing the model. + G (Grid): Holds essential parameters describing the model. """ self.update_magnetic_gpu(np.int32(self.xs), np.int32(self.xf), np.int32(self.ys), np.int32(self.yf), np.int32(self.zs), np.int32(self.zf), np.int32(self.HPhi1_gpu.shape[1]), np.int32(self.HPhi1_gpu.shape[2]), np.int32(self.HPhi1_gpu.shape[3]), np.int32(self.HPhi2_gpu.shape[1]), np.int32(self.HPhi2_gpu.shape[2]), np.int32(self.HPhi2_gpu.shape[3]), np.int32(self.thickness), G.ID_gpu.gpudata, G.Ex_gpu.gpudata, G.Ey_gpu.gpudata, G.Ez_gpu.gpudata, G.Hx_gpu.gpudata, G.Hy_gpu.gpudata, G.Hz_gpu.gpudata, self.HPhi1_gpu.gpudata, self.HPhi2_gpu.gpudata, self.HRA_gpu.gpudata, self.HRB_gpu.gpudata, self.HRE_gpu.gpudata, self.HRF_gpu.gpudata, floattype(self.d), block=G.tpb, grid=self.bpg) def pml_information(G): - # no pml + # No pml if all(value == 0 for value in G.pmlthickness.values()): return 'PML: switched off' @@ -423,7 +419,8 @@ def pml_information(G): for key, value in G.pmlthickness.items(): pmlinfo += '{}: {}, '.format(key, value) pmlinfo = pmlinfo[:-2] + ' cells' - return 'PML: formulation: {}, order: {}, thickness: {}'.format(G.pmlformulation, len(G.cfs), pmlinfo) + return '\nPML: formulation: {}, order: {}, thickness: {}'.format(G.pmlformulation, len(G.cfs), pmlinfo) + def build_pml(G, key, value): """This function builds instances of the PML and calculates the initial @@ -431,56 +428,57 @@ def build_pml(G, key, value): (based on underlying material er and mr from solid array). Args: - G (class): Grid class instance - holds essential parameters - describing the model. + G (Grid): Holds essential parameters describing the model. + key (str): Identifier of PML slab. + value (int): Thickness of PML slab in cells. """ - if value > 0: - sumer = 0 # Sum of relative permittivities in PML slab - summr = 0 # Sum of relative permeabilities in PML slab - if key[0] == 'x': - if key == 'x0': - pml = PML(G, ID=key, direction='xminus', xf=value, yf=G.ny, zf=G.nz) - elif key == 'xmax': - pml = PML(G, ID=key, direction='xplus', xs=G.nx - value, xf=G.nx, yf=G.ny, zf=G.nz) - G.pmls.append(pml) + sumer = 0 # Sum of relative permittivities in PML slab + summr = 0 # Sum of relative permeabilities in PML slab + + if key[0] == 'x': + if key == 'x0': + pml = PML(G, ID=key, direction='xminus', xf=value, yf=G.ny, zf=G.nz) + elif key == 'xmax': + pml = PML(G, ID=key, direction='xplus', xs=G.nx - value, xf=G.nx, yf=G.ny, zf=G.nz) + G.pmls.append(pml) + for j in range(G.ny): + for k in range(G.nz): + numID = G.solid[pml.xs, j, k] + material = next(x for x in G.materials if x.numID == numID) + sumer += material.er + summr += material.mr + averageer = sumer / (G.ny * G.nz) + averagemr = summr / (G.ny * G.nz) + + elif key[0] == 'y': + if key == 'y0': + pml = PML(G, ID=key, direction='yminus', yf=value, xf=G.nx, zf=G.nz) + elif key == 'ymax': + pml = PML(G, ID=key, direction='yplus', ys=G.ny - value, xf=G.nx, yf=G.ny, zf=G.nz) + G.pmls.append(pml) + for i in range(G.nx): + for k in range(G.nz): + numID = G.solid[i, pml.ys, k] + material = next(x for x in G.materials if x.numID == numID) + sumer += material.er + summr += material.mr + averageer = sumer / (G.nx * G.nz) + averagemr = summr / (G.nx * G.nz) + + elif key[0] == 'z': + if key == 'z0': + pml = PML(G, ID=key, direction='zminus', zf=value, xf=G.nx, yf=G.ny) + elif key == 'zmax': + pml = PML(G, ID=key, direction='zplus', zs=G.nz - value, xf=G.nx, yf=G.ny, zf=G.nz) + G.pmls.append(pml) + for i in range(G.nx): for j in range(G.ny): - for k in range(G.nz): - numID = G.solid[pml.xs, j, k] - material = next(x for x in G.materials if x.numID == numID) - sumer += material.er - summr += material.mr - averageer = sumer / (G.ny * G.nz) - averagemr = summr / (G.ny * G.nz) + numID = G.solid[i, j, pml.zs] + material = next(x for x in G.materials if x.numID == numID) + sumer += material.er + summr += material.mr + averageer = sumer / (G.nx * G.ny) + averagemr = summr / (G.nx * G.ny) - elif key[0] == 'y': - if key == 'y0': - pml = PML(G, ID=key, direction='yminus', yf=value, xf=G.nx, zf=G.nz) - elif key == 'ymax': - pml = PML(G, ID=key, direction='yplus', ys=G.ny - value, xf=G.nx, yf=G.ny, zf=G.nz) - G.pmls.append(pml) - for i in range(G.nx): - for k in range(G.nz): - numID = G.solid[i, pml.ys, k] - material = next(x for x in G.materials if x.numID == numID) - sumer += material.er - summr += material.mr - averageer = sumer / (G.nx * G.nz) - averagemr = summr / (G.nx * G.nz) - - elif key[0] == 'z': - if key == 'z0': - pml = PML(G, ID=key, direction='zminus', zf=value, xf=G.nx, yf=G.ny) - elif key == 'zmax': - pml = PML(G, ID=key, direction='zplus', zs=G.nz - value, xf=G.nx, yf=G.ny, zf=G.nz) - G.pmls.append(pml) - for i in range(G.nx): - for j in range(G.ny): - numID = G.solid[i, j, pml.zs] - material = next(x for x in G.materials if x.numID == numID) - sumer += material.er - summr += material.mr - averageer = sumer / (G.nx * G.ny) - averagemr = summr / (G.nx * G.ny) - - pml.calculate_update_coeffs(averageer, averagemr, G) + pml.calculate_update_coeffs(averageer, averagemr, G) diff --git a/gprMax/receivers.py b/gprMax/receivers.py index e7e10bb4..c7ddda43 100644 --- a/gprMax/receivers.py +++ b/gprMax/receivers.py @@ -48,8 +48,7 @@ def gpu_initialise_rx_arrays(G): components for receivers. Args: - G (class): Grid class instance - holds essential parameters - describing the model. + G (Grid): Holds essential parameters describing the model. """ import pycuda.gpuarray as gpuarray @@ -80,8 +79,7 @@ def gpu_get_rx_array(rxs_gpu, rxcoords_gpu, G): rxs_gpu (float): numpy array of receiver data from GPU - rows are field components; columns are iterations; pages are receivers. rxcoords_gpu (float): numpy array of receiver coordinates from GPU. - G (class): Grid class instance - holds essential parameters - describing the model. + G (Grid): Holds essential parameters describing the model. """ for rx in G.rxs: diff --git a/gprMax/scene.py b/gprMax/scene.py index c1833124..3e99d3a7 100644 --- a/gprMax/scene.py +++ b/gprMax/scene.py @@ -1,14 +1,32 @@ -from .user_inputs import create_user_input_points -from .materials import create_built_in_materials +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . + +from .cmds_geometry.cmds_geometry import UserObjectGeometry +from .cmds_geometry.fractal_box_builder import FractalBoxBuilder +from .cmds_multiple import UserObjectMulti from .cmds_single_use import UserObjectSingle from .cmds_single_use import Domain from .cmds_single_use import Discretisation from .cmds_single_use import TimeWindow -from .cmds_multiple import UserObjectMulti -from .subgrids.user_objects import SubGridBase as SubGridUserBase -from .cmds_geometry.cmds_geometry import UserObjectGeometry from .exceptions import CmdInputError -from .cmds_geometry.fractal_box_builder import FractalBoxBuilder +from .materials import create_built_in_materials +from .subgrids.user_objects import SubGridBase as SubGridUserBase +from .user_inputs import create_user_input_points from .utilities import human_size @@ -39,20 +57,20 @@ class Scene: def process_subgrid_commands(self, subgrids): - # check for subgrid user objects + # Check for subgrid user objects def func(obj): if isinstance(obj, SubGridUserBase): return True else: return False - # subgrid user objects + # Subgrid user objects subgrid_cmds = list(filter(func, self.multiple_cmds)) - # iterate through the user command objects under the subgrid user object + # Iterate through the user command objects under the subgrid user object for sg_cmd in subgrid_cmds: - # when the subgrid is created its reference is attached to its user - # object. this reference allows the multi and geo user objects + # When the subgrid is created its reference is attached to its user + # object. This reference allows the multi and geo user objects # to build in the correct subgrid. sg = sg_cmd.subgrid self.process_cmds(sg_cmd.children_multiple, sg) @@ -89,7 +107,6 @@ class Scene: self.process_cmds(cmds_unique, G) - def create_internal_objects(self, G): # fractal box commands have an additional nonuser object which diff --git a/gprMax/solvers.py b/gprMax/solvers.py index 9370a908..44064072 100644 --- a/gprMax/solvers.py +++ b/gprMax/solvers.py @@ -15,17 +15,25 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -from gprMax.updates import CPUUpdates -from gprMax.updates import GPUUpdates -from .subgrids.updates import create_updates as create_subgrid_updates -from gprMax.utilities import timer + +import gprMax.config as config from .grid import FDTDGrid from .grid import GPUGrid -import gprMax.config as config +from .subgrids.updates import create_updates as create_subgrid_updates +from .updates import CPUUpdates +from .updates import GPUUpdates def create_G(sim_config): - """Returns the configured solver.""" + """Create grid object according to solver. + + Args: + sim_config (SimConfig): Simulation level configuration object. + + Returns: + G (Grid): Holds essential parameters describing the model. + """ + if sim_config.gpu: G = GPUGrid() elif sim_config.subgrid: @@ -37,37 +45,46 @@ def create_G(sim_config): def create_solver(G, sim_config): - """Returns the configured solver.""" + """Create configured solver object. + + Args: + G (Grid): Holds essential parameters describing the model. + sim_config (SimConfig): simulation level configuration object. + + Returns: + solver (Solver): solver object. + """ + if sim_config.gpu: updates = GPUUpdates(G) solver = Solver(updates) elif sim_config.subgrid: updates = create_subgrid_updates(G) solver = Solver(updates, hsg=True) + # A large range of different functions exist to advance the time step for + # dispersive materials. The correct function is set here based on the + # the required numerical precision and dispersive material type. + props = updates.adapt_dispersive_config(config) + updates.set_dispersive_updates(props) else: updates = CPUUpdates(G) solver = Solver(updates) - # a large range of function exist to advance the time step for dispersive - # materials. The correct function is set here based on the - # the required numerical precision and dispersive material type. - props = updates.adapt_dispersive_config(config) - updates.set_dispersive_updates(props) + props = updates.adapt_dispersive_config(config) + updates.set_dispersive_updates(props) + return solver - class Solver: - """Generic solver for Update objects""" def __init__(self, updates, hsg=False): - """Context for the model to run in. Sub-class this with contexts - i.e. an MPI context. - - Args: - updates (Updates): updates contains methods to run FDTD algorithm - iterator (iterator): can be range() or tqdm() """ + Args: + updates (Updates): Updates contains methods to run FDTD algorithm. + hsg (bool): Use sub-gridding. + """ + self.updates = updates self.hsg = hsg @@ -75,8 +92,14 @@ class Solver: return self.updates.G def solve(self, iterator): - """Time step the FDTD model.""" - tsolvestart = timer() + """Time step the FDTD model. + + Args: + iterator (iterator): can be range() or tqdm() + """ + + self.updates.time_start() + for iteration in iterator: self.updates.store_outputs() self.updates.store_snapshots(iteration) @@ -92,5 +115,8 @@ class Solver: self.updates.hsg_1() self.updates.update_electric_b() - tsolve = timer() - tsolvestart + self.updates.finalise() + tsolve = self.updates.calculate_tsolve() + self.updates.cleanup() + return tsolve diff --git a/gprMax/sources.py b/gprMax/sources.py index 0bd15189..2067495d 100644 --- a/gprMax/sources.py +++ b/gprMax/sources.py @@ -47,7 +47,7 @@ class Source(object): """Calculates all waveform values for source for duration of simulation. Args: - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ # Waveform values for electric sources - calculated half a timestep later self.waveformvaluesJ = np.zeros((G.iterations), dtype=config.dtypes['float_or_double']) @@ -85,7 +85,7 @@ class VoltageSource(Source): updatecoeffsE (memory view): numpy array of electric field update coefficients. ID (memory view): numpy array of numeric IDs corresponding to materials in the model. Ex, Ey, Ez (memory view): numpy array of electric field values. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: @@ -117,7 +117,7 @@ class VoltageSource(Source): voltage source conductivity to the underlying parameters. Args: - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ if self.resistance != 0: @@ -161,7 +161,7 @@ class HertzianDipole(Source): updatecoeffsE (memory view): numpy array of electric field update coefficients. ID (memory view): numpy array of numeric IDs corresponding to materials in the model. Ex, Ey, Ez (memory view): numpy array of electric field values. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: @@ -193,7 +193,7 @@ class MagneticDipole(Source): updatecoeffsH (memory view): numpy array of magnetic field update coefficients. ID (memory view): numpy array of numeric IDs corresponding to materials in the model. Hx, Hy, Hz (memory view): numpy array of magnetic field values. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: @@ -217,7 +217,7 @@ def gpu_initialise_src_arrays(sources, G): Args: sources (list): List of sources of one class, e.g. HertzianDipoles. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. Returns: srcinfo1_gpu (int): numpy array of source cell coordinates and polarisation information. @@ -260,13 +260,13 @@ def gpu_initialise_src_arrays(sources, G): class TransmissionLine(Source): """A transmission line source is a one-dimensional transmission line - which is attached virtually to a grid cell. + which is attached virtually to a grid cell. """ def __init__(self, G): """ Args: - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ super().__init__() @@ -299,10 +299,11 @@ class TransmissionLine(Source): def calculate_incident_V_I(self, G): """Calculates the incident voltage and current with a long length - transmission line not connected to the main grid from: http://dx.doi.org/10.1002/mop.10415 + transmission line not connected to the main grid + from: http://dx.doi.org/10.1002/mop.10415 Args: - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ for iteration in range(G.iterations): @@ -318,7 +319,7 @@ class TransmissionLine(Source): """Updates absorbing boundary condition at end of the transmission line. Args: - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ h = (config.c * G.dt - self.dl) / (config.c * G.dt + self.dl) @@ -332,7 +333,7 @@ class TransmissionLine(Source): Args: iteration (int): Current iteration (timestep). - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ # Update all the voltage values along the line @@ -349,7 +350,7 @@ class TransmissionLine(Source): Args: iteration (int): Current iteration (timestep). - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ # Update all the current values along the line @@ -366,7 +367,7 @@ class TransmissionLine(Source): updatecoeffsE (memory view): numpy array of electric field update coefficients. ID (memory view): numpy array of numeric IDs corresponding to materials in the model. Ex, Ey, Ez (memory view): numpy array of electric field values. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: @@ -393,7 +394,7 @@ class TransmissionLine(Source): updatecoeffsH (memory view): numpy array of magnetic field update coefficients. ID (memory view): numpy array of numeric IDs corresponding to materials in the model. Hx, Hy, Hz (memory view): numpy array of magnetic field values. - G (class): Grid class instance - holds essential parameters describing the model. + G (Grid): Holds essential parameters describing the model. """ if iteration * G.dt >= self.start and iteration * G.dt <= self.stop: diff --git a/gprMax/subgrids/base.py b/gprMax/subgrids/base.py index 0b1364cd..7a3bd666 100644 --- a/gprMax/subgrids/base.py +++ b/gprMax/subgrids/base.py @@ -15,10 +15,11 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -from ..grid import FDTDGrid -from scipy.constants import c import numpy as np +from scipy.constants import c + +from ..grid import FDTDGrid class SubGridBase(FDTDGrid): diff --git a/gprMax/subgrids/multi.py b/gprMax/subgrids/multi.py index 39b72822..9d92825e 100644 --- a/gprMax/subgrids/multi.py +++ b/gprMax/subgrids/multi.py @@ -15,17 +15,18 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . + from ..receivers import Rx class ReferenceRx(Rx): - """Receiver that micks a receiver in coarse grid.""" - - """We often want to compare an output in a fine reference solution with a - the solution in the coarse grid of a subgridded solution. This receiver - moves the output points in the fine grid such that they are in the same - position as the coarse grid. + """Receiver that mimicks a receiver in coarse grid. + We often want to compare an output in a fine reference solution with a + the solution in the coarse grid of a subgridded solution. This receiver + moves the output points in the fine grid such that they are in the same + position as the coarse grid. """ + def __init__(self): """Constructor.""" super().__init__() @@ -33,19 +34,14 @@ class ReferenceRx(Rx): def get_field(self, str_id, field): """Return the field value at the equivalent coarse yee cell. - Parameters - ---------- - str_id : str - 'Ex' etc... - field : np array - e.g. numpy array of grid.Ez - - Returns - ------- - float - Field value + Args: + str_id (str): 'Ex' etc... + field (nparray): Numpy array of grid.Ez + Returns: + e (float): Field value. """ + d = { 'Ex': self.get_Ex_from_field, 'Ey': self.get_Ey_from_field, @@ -56,21 +52,17 @@ class ReferenceRx(Rx): } e = d[str_id](field) + return e def get_Ex_from_field(self, Ex): """Return the Ex field value from the equivalent coarse yee cell. - Parameters - ---------- - Ex : 3d numpy array - e.g. grid.Ex - - Returns - ------- - float - Ex field value + Args: + Ex (nparray): Numpy array of Ex field. + Returns: + e (float): Ex field value. """ # offset = ratio // 2 @@ -80,16 +72,11 @@ class ReferenceRx(Rx): def get_Ey_from_field(self, Ey): """Return the Ey field value from the equivalent coarse yee cell. - Parameters - ---------- - Ex : 3d numpy array - e.g. grid.Ex - - Returns - ------- - float - Ey field value + Args: + Ey (nparray): Numpy array of Ey field. + Returns: + e (float): Ey field value. """ e = Ey[self.xcoord, self.ycoord + self.offset, self.zcoord] return e @@ -97,33 +84,24 @@ class ReferenceRx(Rx): def get_Ez_from_field(self, Ez): """Return the Ez field value from the equivalent coarse yee cell. - Parameters - ---------- - Ex : 3d numpy array - e.g. grid.Ez - - Returns - ------- - float - Ez field value + Args: + Ez (nparray): Numpy array of Ez field. + Returns: + e (float): Ez field value. """ + e = Ez[self.xcoord, self.ycoord, self.zcoord + self.offset] return e def get_Hx_from_field(self, Hx): """Return the Hx field value from the equivalent coarse yee cell. - Parameters - ---------- - Ex : 3d numpy array - e.g. grid.Hx - - Returns - ------- - float - Hx field value + Args: + Hx (nparray): Numpy array of Hx field. + Returns: + e (float): Hx field value. """ e = Hx[self.xcoord, self.ycoord + self.offset, self.zcoord + self.offset] return e @@ -131,16 +109,11 @@ class ReferenceRx(Rx): def get_Hy_from_field(self, Hy): """Return the Hy field value from the equivalent coarse yee cell. - Parameters - ---------- - Ex : 3d numpy array - e.g. grid.Hx - - Returns - ------- - float - Hy field value + Args: + Hy (nparray): Numpy array of Hy field. + Returns: + e (float): Hy field value. """ e = Hy[self.xcoord + self.offset, self.ycoord, self.zcoord + self.offset] return e @@ -148,16 +121,11 @@ class ReferenceRx(Rx): def get_Hz_from_field(self, Hz): """Return the Hz field value from the equivalent coarse yee cell. - Parameters - ---------- - Ex : 3d numpy array - e.g. grid.Hx - - Returns - ------- - float - Hz field value + Args: + Hz (nparray): Numpy array of Hz field. + Returns: + e (float): Hz field value. """ e = Hz[self.xcoord + self.offset, self.ycoord + self.offset, self.zcoord] return e diff --git a/gprMax/subgrids/precursor_nodes.py b/gprMax/subgrids/precursor_nodes.py index f9d20b66..a14a06b4 100644 --- a/gprMax/subgrids/precursor_nodes.py +++ b/gprMax/subgrids/precursor_nodes.py @@ -15,9 +15,11 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . + +import sys + import numpy as np from scipy import interpolate -import sys def calculate_weighting_coefficients(x1, x): @@ -180,8 +182,7 @@ class PrecusorNodesBase(object): setattr(self, f, val) def calc_exact_field(self, field_names): - """ - Function to set the fields used in update calculations to the + """Function to set the fields used in update calculations to the values at the current main time step. i.e. ey_left = copy.ey_left_1 """ diff --git a/gprMax/subgrids/subgrid_hsg.py b/gprMax/subgrids/subgrid_hsg.py index 47003c96..9c615e48 100644 --- a/gprMax/subgrids/subgrid_hsg.py +++ b/gprMax/subgrids/subgrid_hsg.py @@ -15,15 +15,16 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . + +from colorama import init, Fore, Style +init() + from .base import SubGridBase from ..cython.fields_updates_hsg import cython_update_is from ..cython.fields_updates_hsg import cython_update_magnetic_os from ..cython.fields_updates_hsg import cython_update_electric_os from ..utilities import human_size -from colorama import init, Fore, Style -init() - class SubGridHSG(SubGridBase): @@ -35,9 +36,11 @@ class SubGridHSG(SubGridBase): def update_magnetic_is(self, precursors): """Update the subgrid nodes at the IS with the currents derived - from the main grid. - Args: nwl, nwm, nwn, face, field, inc_field, lookup_id, sign, mod, co + from the main grid. + Args: + nwl, nwm, nwn, face, field, inc_field, lookup_id, sign, mod, co """ + # Hz = c0Hz - c1Ey + c2Ex # Hy = c0Hy - c3Ex + c1Ez # Hx = c0Hx - c2Ez + c3Ey diff --git a/gprMax/subgrids/updates.py b/gprMax/subgrids/updates.py index 99150ef7..77c32e3c 100644 --- a/gprMax/subgrids/updates.py +++ b/gprMax/subgrids/updates.py @@ -15,12 +15,11 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -from ..exceptions import GeneralError -from .subgrid_hsg import SubGridHSG +from ..exceptions import GeneralError from .precursor_nodes import PrecursorNodes from .precursor_nodes import PrecursorNodesFiltered - +from .subgrid_hsg import SubGridHSG from ..updates import CPUUpdates @@ -66,18 +65,17 @@ class SubgridUpdates(CPUUpdates): class SubgridUpdater(CPUUpdates): """Class to handle updating the electric and magnetic fields of an HSG - subgrid. The IS, OS, subgrid region and the electric/magnetic sources are updated - using the precursor regions. + subgrid. The IS, OS, subgrid region and the electric/magnetic sources + are updated using the precursor regions. """ def __init__(self, subgrid, precursors, G): """ Args: - subgrid (SubGrid3d): Subgrid to be updated + subgrid (SubGrid3d): Subgrid to be updated. precursors (PrecursorNodes): Precursor nodes associated with - the subgrid - contain interpolated fields - G (class): Grid class instance - holds essential parameters - describing the model. + the subgrid - contain interpolated fields. + G (Grid): Holds essential parameters describing the model. """ super().__init__(subgrid) self.precursors = precursors @@ -86,7 +84,9 @@ class SubgridUpdater(CPUUpdates): def hsg_1(self): """This is the first half of the subgrid update. Takes the time step - up to the main grid magnetic update""" + up to the main grid magnetic update. + """ + G = self.G sub_grid = self.grid precursors = self.precursors @@ -125,7 +125,9 @@ class SubgridUpdater(CPUUpdates): def hsg_2(self): """This is the first half of the subgrid update. Takes the time step - up to the main grid electric update""" + up to the main grid electric update. + """ + G = self.G sub_grid = self.grid precursors = self.precursors diff --git a/gprMax/subgrids/user_objects.py b/gprMax/subgrids/user_objects.py index 0917bf6d..ca5b7657 100644 --- a/gprMax/subgrids/user_objects.py +++ b/gprMax/subgrids/user_objects.py @@ -15,21 +15,24 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . -from .subgrid_hsg import SubGridHSG as SubGridHSGUser -from .multi import ReferenceRx as ReferenceRxUser -from ..exceptions import CmdInputError -from ..cmds_multiple import UserObjectMulti -from ..cmds_geometry.cmds_geometry import UserObjectGeometry -from ..cmds_multiple import Rx -from gprMax import config from copy import copy import numpy as np +from gprMax import config +from ..cmds_geometry.cmds_geometry import UserObjectGeometry +from ..cmds_multiple import UserObjectMulti +from ..cmds_multiple import Rx +from ..exceptions import CmdInputError +from .multi import ReferenceRx as ReferenceRxUser +from .subgrid_hsg import SubGridHSG as SubGridHSGUser + class SubGridBase(UserObjectMulti): - """Class to allow UserObjectMulti and UserObjectGeometry to be nested in SubGrid type user objects.""" + """Class to allow UserObjectMulti and UserObjectGeometry to be nested + in SubGrid type user objects. + """ def __init__(self, **kwargs): """Constructor.""" @@ -198,7 +201,6 @@ class SubGridHSG(SubGridBase): return sg - class ReferenceRx(Rx): """ReferenceRx User Object.""" diff --git a/gprMax/updates.py b/gprMax/updates.py index 2ab9957b..bf7cb0d7 100644 --- a/gprMax/updates.py +++ b/gprMax/updates.py @@ -15,32 +15,53 @@ # # You should have received a copy of the GNU General Public License # along with gprMax. If not, see . + from importlib import import_module -from gprMax.fields_outputs import store_outputs + import gprMax.config as config -from gprMax.cython.fields_updates_normal import update_electric -from gprMax.cython.fields_updates_normal import update_magnetic +from .cuda.fields_updates import kernel_template_fields +from .cuda.snapshots import kernel_template_store_snapshot +from .cuda.source_updates import kernel_template_sources +from .cython.fields_updates_normal import update_electric +from .cython.fields_updates_normal import update_magnetic +from .fields_outputs import store_outputs +from .receivers import gpu_initialise_rx_arrays +from .receivers import gpu_get_rx_array +from .snapshots import Snapshot +from .snapshots import gpu_initialise_snapshot_array +from .snapshots import gpu_get_snapshot_array +from .sources import gpu_initialise_src_arrays +from .utilities import timer class CPUUpdates: + """Defines update functions for CPU-based solver.""" def __init__(self, G): + """ + Args: + G (FDTDGrid): FDTD grid object + """ self.grid = G self.dispersive_update_a = None self.dispersive_update_b = None def store_outputs(self): - # Store field component values for every receiver and transmission line + """Store field component values for every receiver and transmission line.""" store_outputs(self.grid) def store_snapshots(self, iteration): - # Store any snapshots + """Store any snapshots. + + Args: + iteration (int): iteration number. + """ for snap in self.grid.snapshots: if snap.time == iteration + 1: snap.store(self.grid) def update_magnetic(self): - # Update magnetic field components + """Update magnetic field components.""" update_magnetic(self.grid.nx, self.grid.ny, self.grid.nz, @@ -55,12 +76,12 @@ class CPUUpdates: self.grid.Hz) def update_magnetic_pml(self): - # Update magnetic field components with the PML correction + """Update magnetic field components with the PML correction.""" for pml in self.grid.pmls: pml.update_magnetic(self.grid) def update_magnetic_sources(self): - # Update magnetic field components from sources + """Update magnetic field components from sources.""" for source in self.grid.transmissionlines + self.grid.magneticdipoles: source.update_magnetic(self.grid.iteration, self.grid.updatecoeffsH, @@ -71,8 +92,8 @@ class CPUUpdates: self.grid) def update_electric_a(self): - # Update electric field components - # All materials are non-dispersive so do standard update + """Update electric field components.""" + # All materials are non-dispersive so do standard update. if config.materials['maxpoles'] == 0: update_electric(self.grid.nx, self.grid.ny, @@ -109,21 +130,24 @@ class CPUUpdates: self.grid.Hz) def update_electric_pml(self): - # Update electric field components with the PML correction + """Update electric field components with the PML correction.""" for pml in self.grid.pmls: pml.update_electric(self.grid) def update_electric_sources(self): - # Update electric field components from sources (update any Hertzian dipole sources last) + """Update electric field components from sources - + update any Hertzian dipole sources last. + """ for source in self.grid.voltagesources + self.grid.transmissionlines + self.grid.hertziandipoles: source.update_electric(self.grid.iteration, self.grid.updatecoeffsE, self.grid.ID, self.grid.Ex, self.grid.Ey, self.grid.Ez, self.grid) self.grid.iteration += 1 def update_electric_b(self): - # If there are any dispersive materials do 2nd part of dispersive update - # (it is split into two parts as it requires present and updated electric - # field values). Therefore it can only be completely updated after the - # electric field has been updated by the PML and source updates. + """If there are any dispersive materials do 2nd part of dispersive update - + it is split into two parts as it requires present and updated electric + field values. Therefore it can only be completely updated after the + electric field has been updated by the PML and source updates. + """ if config.materials['maxpoles'] != 0: self.dispersive_update_b(self.grid.nx, self.grid.ny, @@ -140,16 +164,21 @@ class CPUUpdates: self.grid.Ez) def adapt_dispersive_config(self, config): + """Set properties for disperive materials. + Args: + config (): + + Returns: + props (Props): Dispersive material properties. + """ if config.materials['maxpoles'] > 1: poles = 'multi' - else: poles = '1' if config.precision == 'single': type = 'float' - else: type = 'double' @@ -169,8 +198,11 @@ class CPUUpdates: return props def set_dispersive_updates(self, props): - """Function to set dispersive update functions based on model.""" + """Set dispersive update functions. + Args: + props (Props): Dispersive material properties. + """ update_f = 'update_electric_dispersive_{}pole_{}_{}_{}' disp_a = update_f.format(props.poles, 'A', props.precision, props.dispersion_type) disp_b = update_f.format(props.poles, 'B', props.precision, props.dispersion_type) @@ -181,6 +213,367 @@ class CPUUpdates: self.dispersive_update_a = disp_a_f self.dispersive_update_b = disp_b_f + def time_start(self): + self.timestart = timer() + + def calculate_tsolve(self): + return timer() - self.timestart + + def finalise(self): + pass + + def cleanup(self): + pass + class GPUUpdates: - pass + """Defines update functions for GPU-based (CUDA) solver.""" + + def __init__(self, G): + """ + Args: + G (FDTDGrid): FDTD grid object + """ + + import pycuda.driver as drv + from pycuda.compiler import SourceModule + drv.init() + + # Suppress nvcc warnings on Windows + if sys.platform == 'win32': + self.compiler_opts = ['-w'] + else: + self.compiler_opts = None + + # Create device handle and context on specifc GPU device (and make it current context) + self.dev = drv.Device(self.grid.gpu.deviceID) + self.ctx = dev.make_context() + + self.grid = G + self.dispersive_update_a = None + self.dispersive_update_b = None + + # Initialise arrays on GPU, prepare kernels, and get kernel functions + self.set_field_kernels() + self.set_pml_kernels() + self.set_rx_kernel() + self.set_src_kernels() + self.set_snapshot_kernel() + + def set_field_kernels(self): + # Electric and magnetic field updates - prepare kernels, and get kernel functions + if config.materials['maxpoles'] > 0: + kernels_fields = SourceModule(kernels_template_fields.substitute(REAL=cudafloattype, COMPLEX=cudacomplextype, N_updatecoeffsE=self.grid.updatecoeffsE.size, N_updatecoeffsH=self.grid.updatecoeffsH.size, NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], NY_MATDISPCOEFFS=self.grid.updatecoeffsdispersive.shape[1], NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1, NX_ID=self.grid.ID.shape[1], NY_ID=self.grid.ID.shape[2], NZ_ID=self.grid.ID.shape[3], NX_T=self.grid.Tx.shape[1], NY_T=self.grid.Tx.shape[2], NZ_T=self.grid.Tx.shape[3]), options=self.compiler_opts) + else: # Set to one any substitutions for dispersive materials + kernels_fields = SourceModule(kernels_template_fields.substitute(REAL=cudafloattype, COMPLEX=cudacomplextype, N_updatecoeffsE=self.grid.updatecoeffsE.size, N_updatecoeffsH=self.grid.updatecoeffsH.size, NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], NY_MATDISPCOEFFS=1, NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1, NX_ID=self.grid.ID.shape[1], NY_ID=self.grid.ID.shape[2], NZ_ID=self.grid.ID.shape[3], NX_T=1, NY_T=1, NZ_T=1), options=self.compiler_opts) + self.update_electric = kernels_fields.get_function("update_electric") + self.update_magnetic = kernels_fields.get_function("update_magnetic") + if self.grid.updatecoeffsE.nbytes + self.grid.updatecoeffsH.nbytes > self.grid.gpu.constmem: + raise GeneralError('Too many materials in the model to fit onto constant memory of size {} on {} - {} GPU'.format(human_size(self.grid.gpu.constmem), self.grid.gpu.deviceID, self.grid.gpu.name)) + self.copy_mat_coeffs() + + # Electric and magnetic field updates - dispersive materials - get kernel functions and initialise array on GPU + if config.materials['maxpoles'] > 0: # If there are any dispersive materials (updates are split into two parts as they require present and updated electric field values). + self.dispersive_update_a = kernels_fields.get_function("update_electric_dispersive_A") + self.dispersive_update_b = kernels_fields.get_function("update_electric_dispersive_B") + self.grid.gpu_initialise_dispersive_arrays() + + # Electric and magnetic field updates - set blocks per grid and initialise field arrays on GPU + self.grid.gpu_set_blocks_per_grid() + self.grid.gpu_initialise_arrays() + + def set_pml_kernels(self): + if self.grid.pmls: + # PMLS - prepare kernels and get kernel functions + pmlmodulelectric = 'gprMax.cuda.pml_updates_electric_' + self.grid.pmlformulation + kernelelectricfunc = getattr(import_module(pmlmodulelectric), 'kernels_template_pml_electric_' + self.grid.pmlformulation) + pmlmodulemagnetic = 'gprMax.cuda.pml_updates_magnetic_' + self.grid.pmlformulation + kernelmagneticfunc = getattr(import_module(pmlmodulemagnetic), 'kernels_template_pml_magnetic_' + self.grid.pmlformulation) + kernels_pml_electric = SourceModule(kernelelectricfunc.substitute(REAL=cudafloattype, N_updatecoeffsE=self.grid.updatecoeffsE.size, NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1, NX_ID=self.grid.ID.shape[1], NY_ID=self.grid.ID.shape[2], NZ_ID=self.grid.ID.shape[3]), options=self.compiler_opts) + kernels_pml_magnetic = SourceModule(kernelmagneticfunc.substitute(REAL=cudafloattype, N_updatecoeffsH=self.grid.updatecoeffsH.size, NY_MATCOEFFS=self.grid.updatecoeffsH.shape[1], NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1, NX_ID=self.gridG.ID.shape[1], NY_ID=self.grid.ID.shape[2], NZ_ID=self.grid.ID.shape[3]), options=self.compiler_opts) + self.copy_mat_coeffs() + # Set block per grid, initialise arrays on GPU, and get kernel functions + for pml in self.grid.pmls: + pml.gpu_initialise_arrays() + pml.gpu_get_update_funcs(kernels_pml_electric, kernels_pml_magnetic) + pml.gpu_set_blocks_per_grid(self.grid) + + def set_rx_kernel(self): + # Receivers - initialise arrays on GPU, prepare kernel and get kernel function + if self.grid.rxs: + rxcoords_gpu, rxs_gpu = gpu_initialise_rx_arrays(self.grid) + kernel_store_outputs = SourceModule(kernel_template_store_outputs.substitute(REAL=cudafloattype, NY_RXCOORDS=3, NX_RXS=6, NY_RXS=self.grid.iterations, NZ_RXS=len(self.grid.rxs), NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1), options=self.compiler_opts) + self.store_outputs = kernel_store_outputs.get_function("store_outputs") + + def set_src_kernels(self): + # Sources - initialise arrays on GPU, prepare kernel and get kernel function + if self.grid.voltagesources + self.grid.hertziandipoles + self.grid.magneticdipoles: + kernels_sources = SourceModule(kernels_template_sources.substitute(REAL=cudafloattype, N_updatecoeffsE=self.grid.updatecoeffsE.size, N_updatecoeffsH=self.grid.updatecoeffsH.size, NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1], NY_SRCINFO=4, NY_SRCWAVES=self.grid.iterations, NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1, NX_ID=self.grid.ID.shape[1], NY_ID=self.grid.ID.shape[2], NZ_ID=self.grid.ID.shape[3]), options=self.compiler_opts) + self.copy_mat_coeffs() + if self.grid.hertziandipoles: + self.srcinfo1_hertzian_gpu, self.srcinfo2_hertzian_gpu, self.srcwaves_hertzian_gpu = gpu_initialise_src_arrays(self.grid.hertziandipoles, self.grid) + self.update_hertzian_dipole_gpu = kernels_sources.get_function("update_hertzian_dipole") + if self.grid.magneticdipoles: + self.srcinfo1_magnetic_gpu, self.srcinfo2_magnetic_gpu, self.srcwaves_magnetic_gpu = gpu_initialise_src_arrays(self.grid.magneticdipoles, self.grid) + self.update_magnetic_dipole_gpu = kernels_sources.get_function("update_magnetic_dipole") + if self.grid.voltagesources: + self.srcinfo1_voltage_gpu, self.srcinfo2_voltage_gpu, self.srcwaves_voltage_gpu = gpu_initialise_src_arrays(self.grid.voltagesources, self.grid) + self.update_voltage_source_gpu = kernels_sources.get_function("update_voltage_source") + + def set_snapshot_kernel(self): + if self.grid.snapshots: + # Snapshots - initialise arrays on GPU, prepare kernel and get kernel function + self.snapEx_gpu, self.snapEy_gpu, self.snapEz_gpu, self.snapHx_gpu, self.snapHy_gpu, self.snapHz_gpu = gpu_initialise_snapshot_array(self.grid) + kernel_store_snapshot = SourceModule(kernel_template_store_snapshot.substitute(REAL=cudafloattype, NX_SNAPS=Snapshot.nx_max, NY_SNAPS=Snapshot.ny_max, NZ_SNAPS=Snapshot.nz_max, NX_FIELDS=self.grid.nx + 1, NY_FIELDS=self.grid.ny + 1, NZ_FIELDS=self.grid.nz + 1), options=self.compiler_opts) + self.store_snapshot_gpu = kernel_store_snapshot.get_function("store_snapshot") + + def copy_mat_coeffs(self): + # Copy material coefficient arrays to constant memory of GPU (must be <64KB) + updatecoeffsE = kernels_sources.get_global('updatecoeffsE')[0] + updatecoeffsH = kernels_sources.get_global('updatecoeffsH')[0] + self.drv.memcpy_htod(updatecoeffsE, self.grid.updatecoeffsE) + self.drv.memcpy_htod(updatecoeffsH, self.grid.updatecoeffsH) + + def store_outputs(self): + """Store field component values for every receiver and transmission line.""" + if self.grid.rxs: + self.store_outputs(np.int32(len(self.grid.rxs)), + np.int32(self.grid.iteration), + self.rxcoords_gpu.gpudata, + self.rxs_gpu.gpudata, + self.grid.Ex_gpu.gpudata, + self.grid.Ey_gpu.gpudata, + self.grid.Ez_gpu.gpudata, + self.grid.Hx_gpu.gpudata, + self.grid.Hy_gpu.gpudata, + self.grid.Hz_gpu.gpudata, + block=(1, 1, 1), + grid=(round32(len(self.grid.rxs)), 1, 1)) + + def store_snapshots(self, iteration): + """Store any snapshots. + + Args: + iteration (int): iteration number. + """ + + for i, snap in enumerate(self.grid.snapshots): + if snap.time == iteration + 1: + snapno = 0 if self.grid.snapsgpu2cpu else i + self.store_snapshot_gpu(np.int32(snapno), + np.int32(snap.xs), + np.int32(snap.xf), + np.int32(snap.ys), + np.int32(snap.yf), + np.int32(snap.zs), + np.int32(snap.zf), + np.int32(snap.dx), + np.int32(snap.dy), + np.int32(snap.dz), + self.grid.Ex_gpu.gpudata, + self.grid.Ey_gpu.gpudata, + self.grid.Ez_gpu.gpudata, + self.grid.Hx_gpu.gpudata, + self.grid.Hy_gpu.gpudata, + self.grid.Hz_gpu.gpudata, + self.snapEx_gpu.gpudata, + self.snapEy_gpu.gpudata, + self.snapEz_gpu.gpudata, + self.snapHx_gpu.gpudata, + self.snapHy_gpu.gpudata, + self.snapHz_gpu.gpudata, + block=Snapshot.tpb, + grid=Snapshot.bpg) + if self.grid.snapsgpu2cpu: + gpu_get_snapshot_array(self.grid.snapEx_gpu.get(), + self.grid.snapEy_gpu.get(), + self.grid.snapEz_gpu.get(), + self.grid.snapHx_gpu.get(), + self.grid.snapHy_gpu.get(), + self.grid.snapHz_gpu.get(), + 0, + snap) + + def update_magnetic(self): + """Update magnetic field components.""" + + self.update_magnetic(np.int32(self.grid.nx), + np.int32(self.grid.ny), + np.int32(self.grid.nz), + self.grid.ID_gpu, + self.grid.Hx_gpu, + self.grid.Hy_gpu, + self.grid.Hz_gpu, + self.grid.Ex_gpu, + self.grid.Ey_gpu, + self.grid.Ez_gpu, + block=self.grid.tpb, + grid=self.grid.bpg) + + def update_magnetic_pml(self): + """Update magnetic field components with the PML correction.""" + + for pml in self.grid.pmls: + pml.gpu_update_magnetic(self.grid) + + def update_magnetic_sources(self): + """Update magnetic field components from sources.""" + + if self.grid.magneticdipoles: + self.update_magnetic_dipole_gpu(np.int32(len(self.grid.magneticdipoles)), + np.int32(self.grid.iteration), + config.dtypes['float_or_double'](self.grid.dx), + config.dtypes['float_or_double'](self.grid.dy), + config.dtypes['float_or_double'](self.grid.dz), + self.srcinfo1_magnetic_gpu.gpudata, + self.srcinfo2_magnetic_gpu.gpudata, + self.srcwaves_magnetic_gpu.gpudata, + self.grid.ID_gpu, + self.grid.Hx_gpu, + self.grid.Hy_gpu, + self.grid.Hz_gpu, + block=(1, 1, 1), + grid=(round32(len(self.grid.magneticdipoles)), 1, 1)) + + def update_electric_a(self): + """Update electric field components.""" + + # All materials are non-dispersive so do standard update. + if config.materials['maxpoles'] == 0: + self.update_electric(np.int32(self.grid.nx), + np.int32(self.grid.ny), + np.int32(self.grid.nz), + self.grid.ID_gpu, + self.grid.Ex_gpu, + self.grid.Ey_gpu, + self.grid.Ez_gpu, + self.grid.Hx_gpu, + self.grid.Hy_gpu, + self.grid.Hz_gpu, + block=self.grid.tpb, + grid=self.grid.bpg) + + # If there are any dispersive materials do 1st part of dispersive update + # (it is split into two parts as it requires present and updated electric field values). + else: + self.dispersive_update_a(np.int32(self.grid.nx), + np.int32(self.grid.ny), + np.int32(self.grid.nz), + np.int32(config.materials['maxpoles']), + self.grid.updatecoeffsdispersive_gpu.gpudata, + self.grid.Tx_gpu, + self.grid.Ty_gpu, + self.grid.Tz_gpu, + self.grid.ID_gpu, + self.grid.Ex_gpu, + self.grid.Ey_gpu, + self.grid.Ez_gpu, + self.grid.Hx_gpu, + self.grid.Hy_gpu, + self.grid.Hz_gpu, + block=self.grid.tpb, + grid=self.grid.bpg) + + def update_electric_pml(self): + """Update electric field components with the PML correction.""" + + for pml in self.grid.pmls: + pml.gpu_update_electric(self.grid) + + def update_electric_sources(self): + """Update electric field components from sources - + update any Hertzian dipole sources last. + """ + + if self.grid.voltagesources: + self.update_voltage_source_gpu(np.int32(len(self.grid.voltagesources)), + np.int32(self.grid.iteration), + config.dtypes['float_or_double'](self.grid.dx), + config.dtypes['float_or_double'](self.grid.dy), + config.dtypes['float_or_double'](self.grid.dz), + self.srcinfo1_voltage_gpu.gpudata, + self.srcinfo2_voltage_gpu.gpudata, + self.srcwaves_voltage_gpu.gpudata, + self.grid.ID_gpu, + self.grid.Ex_gpu, + self.grid.Ey_gpu, + self.grid.Ez_gpu, + block=(1, 1, 1), + grid=(round32(len(self.grid.voltagesources)), 1, 1)) + + if self.grid.hertziandipoles: + self.update_hertzian_dipole_gpu(np.int32(len(self.grid.hertziandipoles)), + np.int32(self.grid.iteration), + config.dtypes['float_or_double'](self.grid.dx), + config.dtypes['float_or_double'](self.grid.dy), + config.dtypes['float_or_double'](self.grid.dz), + self.srcinfo1_hertzian_gpu.gpudata, + self.srcinfo2_hertzian_gpu.gpudata, + self.srcwaves_hertzian_gpu.gpudata, + self.grid.ID_gpu, + self.grid.Ex_gpu, + self.grid.Ey_gpu, + self.grid.Ez_gpu, + block=(1, 1, 1), + grid=(round32(len(self.grid.hertziandipoles)), 1, 1)) + + self.grid.iteration += 1 + + def update_electric_b(self): + """If there are any dispersive materials do 2nd part of dispersive update - + it is split into two parts as it requires present and updated electric + field values. Therefore it can only be completely updated after the + electric field has been updated by the PML and source updates. + """ + + if config.materials['maxpoles'] != 0: + self.dispersive_update_b(np.int32(self.grid.nx), + np.int32(self.grid.ny), + np.int32(self.grid.nz), + np.int32(config.materials['maxpoles']), + self.grid.updatecoeffsdispersive_gpu.gpudata, + self.grid.Tx_gpu, + self.grid.Ty_gpu, + self.grid.Tz_gpu, + self.grid.ID_gpu, + self.grid.Ex_gpu, + self.grid.Ey_gpu, + self.grid.Ez_gpu, + block=self.grid.tpb, + grid=self.grid.bpg) + + def time_start(self): + self.iterstart = self.drv.Event() + self.iterend = self.drv.Event() + self.iterstart.record() + + def calculate_tsolve(self): + self.iterend.record() + self.iterend.synchronize() + tsolve = self.iterstart.time_till(self.iterend) * 1e-3 + + return tsolve + + def finalise(self): + # Copy output from receivers array back to correct receiver objects + if self.grid.rxs: + gpu_get_rx_array(self.rxs_gpu.get(), + self.rxcoords_gpu.get(), + self.grid) + + # Copy data from any snapshots back to correct snapshot objects + if self.grid.snapshots and not self.grid.snapsgpu2cpu: + for i, snap in enumerate(self.grid.snapshots): + gpu_get_snapshot_array(self.snapEx_gpu.get(), + self.snapEy_gpu.get(), + self.snapEz_gpu.get(), + self.snapHx_gpu.get(), + self.snapHy_gpu.get(), + self.snapHz_gpu.get(), + i, + snap) + + def cleanup(self): + # Remove context from top of stack and delete + self.ctx.pop() + del self.ctx diff --git a/gprMax/user_inputs.py b/gprMax/user_inputs.py index f25c1e63..67a93fd9 100644 --- a/gprMax/user_inputs.py +++ b/gprMax/user_inputs.py @@ -1,19 +1,37 @@ -from .exceptions import CmdInputError -from .subgrids.base import SubGridBase -from .utilities import round_value -import gprMax.config as config +# Copyright (C) 2015-2019: The University of Edinburgh +# Authors: Craig Warren and Antonis Giannopoulos +# +# This file is part of gprMax. +# +# gprMax is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# gprMax is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with gprMax. If not, see . -import numpy as np from colorama import init from colorama import Fore from colorama import Style - init() +import numpy as np + +import gprMax.config as config +from .exceptions import CmdInputError +from .subgrids.base import SubGridBase +from .utilities import round_value + """Module contains classes to handle points supplied by a user. The classes implement a common interface such that geometry building objects -such as Box or triangle do not need to have any knowledge which grid to which they -are rounding continuous points or checking the point is within the grid. +such as box or triangle do not need to have any knowledge which grid to which +they are rounding continuous points or checking the point is within the grid. Additionally all logic related to rounding points etc is encapulsated here. """ @@ -36,7 +54,6 @@ def create_user_input_points(grid, user_obj): class UserInput: - """Base class to handle (x, y, z) points supplied by the user.""" def __init__(self, grid): @@ -113,8 +130,9 @@ class MainGridUserInput(UserInput): class SubgridUserInput(MainGridUserInput): """Class to handle (x, y, z) points supplied by the user in the sub grid. - This class autotranslates points from main grid to subgrid equivalent (within IS). Useful - if material traverse is not required.""" + This class autotranslates points from main grid to subgrid equivalent (within IS). + Useful if material traverse is not required. + """ def __init__(self, grid): super().__init__(grid) diff --git a/gprMax/utilities.py b/gprMax/utilities.py index 2bd597ef..71630ee0 100644 --- a/gprMax/utilities.py +++ b/gprMax/utilities.py @@ -420,8 +420,10 @@ def timer(): """Function to return the current process wide time in fractional seconds.""" return perf_counter() -class Printer(): +class Printer(): + """Printing information messages.""" + def __init__(self, config): self.printing = config.is_messages() diff --git a/temp/gprMax b/temp/gprMax deleted file mode 160000 index 411a3f9e..00000000 --- a/temp/gprMax +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 411a3f9ec8da465271cf8e75057a597fe52369d4