publicImage/gprMax
1
0
派生 0
镜像自地址 https://gitee.com/sunhf/gprMax.git 已同步 2025-08-07 15:10:13 +08:00
代码 工单 软件包 项目 版本发布 百科 动态

Removal of Taguchi optimisation. Move towards high-level config module, rather than settings in Grid class.

浏览代码
这个提交包含在：
Craig Warren
2019-05-24 09:54:32 +01:00
父节点 142b9b4064
当前提交 61e9d43e26
共有 35 个文件被更改，包括 349 次插入 和 1525 次删除
下载 Patch 文件 下载 Diff 文件 展开所有文件 折叠所有文件

1
README.rst
查看文件

@@ -183,7 +183,6 @@ Argument name Type Description
``-benchmark`` flag 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 <http://docs.gprmax.com/en/latest/benchmarking.html>`_
``--geometry-only`` flag build a model and produce any geometry views but do not run the simulation, e.g. to check the geometry of a model is correct: ``(gprMax)$ python -m gprMax user_models/heterogeneous_soil.in --geometry-only``
``--geometry-fixed`` flag run a series of models where the geometry does not change between models, e.g. a B-scan where *only* the position of simple sources and receivers, moved using ``#src_steps`` and ``#rx_steps``, changes between models.
``--opt-taguchi`` flag run a series of models using an optimisation process based on Taguchi's method. For further details see the `user libraries section of the User Guide <http://docs.gprmax.com/en/latest/user_libs_opt_taguchi.html>`_
``--write-processed`` flag write another input file after any Python code and include commands in the original input file have been processed. Useful for checking that any Python code is being correctly processed into gprMax commands.
``-h`` or ``--help`` flag used to get help on command line options.
====================== ========= ===========

1
docs/source/index.rst
查看文件

@@ -41,7 +41,6 @@ gprMax User Guide
user_libs_antenna_patterns
user_libs_austinman
user_libs_materials
user_libs_opt_taguchi
.. toctree::
:maxdepth: 2

149
docs/source/user_libs_opt_taguchi.rst
查看文件

@@ -1,149 +0,0 @@
User libraries is a sub-package where useful Python modules contributed by users are stored.
*******************************
Optimisation - Taguchi's method
*******************************
Information
===========
**Author/Contact**: Craig Warren (craig.warren@northumbria.ac.uk), Northumbria University, UK
**License**: `Creative Commons Attribution-ShareAlike 4.0 International License <http://creativecommons.org/licenses/by-sa/4.0/>`_
**Attribution/cite**: Warren, C., Giannopoulos, A. (2011). Creating finite-difference time-domain models of commercial ground-penetrating radar antennas using Taguchi's optimization method. *Geophysics*, 76(2), G37-G47. (http://dx.doi.org/10.1190/1.3548506)
The package features an optimisation technique based on Taguchi's method. It allows users to define parameters in an input file and optimise their values based on a fitness function, for example it can be used to optimise material properties or geometry in a simulation.
.. warning::
This package combines a number of advanced features and should not be used without knowledge and familiarity of the underlying techniques. It requires:
* Knowledge of Python to contruct an input file to use with the optimisation
* Familiarity of optimisation techniques, and in particular Taguchi's method
* Careful sanity checks to be made throughout the process
Taguchi's method
----------------
Taguchi's method is based on the concept of the Orthogonal Array (OA) and has the following advantages:
* Simple to implement
* Effective in reduction of experiments
* Fast convergence speed
* Global optimum results
* Independence from initial values of optimisation parameters
Details of Taguchi's method in the context of electromagnetics can be found in [WEN2007a]_ and [WEN2007b]_.
Package overview
================
.. code-block:: none
antenna_bowtie_opt.in
fitness_functions.py
OA_9_4_3_2.npy
OA_18_7_3_2.npy
plot_results.py
* ``antenna_bowtie_opt.in`` is a example model of a bowtie antenna where values of loading resistors are optimised.
* ``fitness_functions.py`` is a module containing fitness functions. There are some pre-built ones but users should add their own here.
* ``OA_9_4_3_2.npy`` and ``OA_18_7_3_2.npy`` are NumPy archives `containing pre-built OAs <http://neilsloane.com/oadir/>`_
* ``plot_results.py`` is a module for plotting the results, such as parameter values and convergence history, from an optimisation process when it has completed.
Implementation
--------------
The process by which Taguchi's method optimises parameters is illustrated in the following figure:
.. figure:: images/user_libs/taguchi_process.png
:width: 300 px
Process associated with Taguchi's method.
In stage 1a, one of the 2 pre-built OAs will automatically be chosen depending on the number of parameters to optimise. Currently, up to 7 independent parameters can be optimised, although a method to construct OAs of any size is under testing.
In stage 1b, a fitness function is required to set a goal against which to compare results from the optimisation process. A number of pre-built fitness functions can be found in the ``fitness_functions.py`` module, e.g. ``minvalue``, ``maxvalue`` and ``xcorr``. Users can also easily add their own fitness functions to this module. All fitness functions must take two arguments:
* ``filename`` a string containing the full path and filename of the output file
* ``args`` a dictionary which can contain any number of additional arguments for the function, e.g. names of outputs (rxs) in the model
Additionally, all fitness functions must return a single fitness value which the optimsation process will aim to maximise.
Stages 2-6 are iterated through by the optimisation process.
Parameters and settings for the optimisation process are specified within a special Python block defined by ``#taguchi`` and ``#end_taguchi`` in the input file. The parameters to optimise must be defined in a dictionary named ``optparams`` and their initial ranges specified as lists with lower and upper values. The fitness function, it's parameters, and a stopping value are defined in dictionary named ``fitness`` which has keys for:
* ``name`` a string that is the name of the fitness function to be used
* ``args`` a dictionary containing arguments to be passed to the fitness function. Within ``args`` there must be a key called ``outputs`` which contains a string or list of the names of one or more outputs (rxs) in the model
* ``stop`` a value from the fitness function which when exceeded the optimisation should stop
Optionally a variable called ``maxiterations`` maybe specified which will set a maximum number of iterations after which the optimisation process will terminate irrespective of any other criteria. If it is not specified it defaults to a maximum of 20 iterations.
There is also a builtin criterion to terminate the optimisation process is successive fitness values are within 0.1% of one another.
How to use the package
======================
The package requires ``#python`` and ``#end_python`` to be used in the input file, as well as ``#taguchi`` and ``#end_taguchi`` for specifying parameters and setting for the optimisation process. A Taguchi optimisation is run using the command line option ``--opt-taguchi``.
Example
-------
The following example demonstrates using the Taguchi optimisation process to optimise values of loading resistors used in a bowtie antenna. The example is slighty contrived as the goal is simply to find values for the resistors that produces a maximum absolute amplitude in the response from the antenna. We already know this should occur when the values of the resistors are at a minimum. Nevertheless, it is useful to illustrate the optimisation process and how to use it.
.. figure:: images/user_libs/antenna_bowtie_opt.png
:width: 600 px
FDTD geometry mesh showing bowtie antenna with slots and loading resistors.
The bowtie design features three vertical slots (y-direction) in each arm of the bowtie. Each slot has different loading resistors, but within each slot there are four resistors of the same value. A resistor is modelled as two parallel edges of a cell. The bowtie is placed on a lossless substrate of relative permittivity 4.8. The antenna is modelled in free space, and an output point of the electric field (named ``Ex60mm``) is specified at a distance of 60mm from the feed of the bowtie (red coloured cell).
.. literalinclude:: ../../user_libs/optimisation_taguchi/antenna_bowtie_opt.in
:language: none
:linenos:
The first part of the input file (lines 1-6) contains the parameters to optimise, their initial ranges, and fitness function information for the optimisation process. Three parameters representing the resistor values are defined with ranges between 0.1 :math:`\Omega` and 1 :math:`k\Omega`. A pre-built fitness function called ``min_max_value`` is specified with a stopping criterion of 10V/m. Arguments for the ``min_max_value`` function are ``type`` given as ``absmax``, i.e. the maximum absolute values, and the output point in the model that will be used in the optimisation is specified as having the name ``Ex60mm``.
The next part of the input file (lines 8-92) contains the model. For the most part there is nothing special about the way the model is defined - a mixture of Python, NumPy and functional forms of the input commands (available by importing the module ``input_cmd_funcs``) are used. However, it is worth pointing out how the values of the parameters to optimise are accessed. On line 29 a NumPy array of the values of the resistors is created. The values are accessed using their names as keys to the ``optparams`` dictionary. On line 30 the values of the resistors are converted to conductivities, which are used to create new materials (line 34-35). The resistors are then built by applying the materials to cell edges (e.g. lines 55-62). The output point in the model in specifed with the name ``Ex60mm`` and as having only an ``Ex`` field output (line 42).
The optimisation process is run on the model using the ``--opt-taguchi`` command line flag.
.. code-block:: none
python -m gprMax user_libs/optimisation_taguchi/antenna_bowtie_opt.in --opt-taguchi
Results
^^^^^^^
When the optimisation has completed a summary will be printed showing histories of the parameter values and the fitness metric. These values are also saved (pickled) to file and can be plotted using the ``plot_results.py`` module, for example:
.. code-block:: none
python -m user_libs.optimisation_taguchi.plot_results user_libs/optimisation_taguchi/antenna_bowtie_opt_hist.pickle
.. code-block:: none
Optimisations summary for: antenna_bowtie_opt_hist.pickle
Number of iterations: 4
History of fitness values: [4.2720928, 5.68856, 5.7023263, 5.7023263]
History of parameter values:
resinner [250.07498, 0.87031555, 0.1, 0.1]
resmiddle [250.07498, 0.87031555, 0.1, 0.1]
resouter [250.07498, 0.87031555, 0.1, 0.1]
.. figure:: images/user_libs/taguchi_fitness_hist.png
:width: 600 px
History of values of fitness metric (absolute maximum).
.. figure:: images/user_libs/taguchi_parameter_hist.png
:width: 600 px
History of values of parameters - ``resinner``, ``resmiddle``, and ``resouter`` (in this case they are all identical).
The optimisation process terminated after 4 iterations because succcessive fitness values were within 0.1% of one another. A maximum absolute amplitude value of 5.7 V/m was achieved when the three resistors had values of 0.1 :math:`\Omega`.

17
gprMax/constants.pxd → gprMax/config.pxd
查看文件

@@ -19,16 +19,11 @@
import numpy as np
cimport numpy as np
# Data types:
# Solid and ID arrays use 32-bit integers (0 to 4294967295)
# Rigid arrays use 8-bit integers (the smallest available type to store true/false)
# Fractal and dispersive coefficient arrays use complex numbers (complextype) which are represented as two floats
# Main field arrays use floats (floattype) and complex numbers (complextype)
# Single precision
ctypedef np.float32_t floattype_t
ctypedef np.complex64_t complextype_t
ctypedef fused floattype_t:
np.float32_t
np.float64_t
# Double precision
# ctypedef np.float64_t floattype_t
# ctypedef np.complex128_t complextype_t
ctypedef fused complextype_t:
np.complex64_t
np.complex128_t

78
gprMax/config.py 普通文件
查看文件

@@ -0,0 +1,78 @@
# 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 <http://www.gnu.org/licenses/>.
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 gprMax.utilities import get_host_info
# Impedance of free space (Ohms)
z0 = np.sqrt(m0 / e0)
# Setting whether messages and progress bars are printed
messages = True
progressbars = messages
# Store information about host machine
hostinfo = get_host_info()
# Store information about any GPUs as a list of GPU objects
gpus = None
# Copy snapshot data from GPU to CPU during simulation
# N.B. This will happen if the requested snapshots are too large to fit
# on the memory of the GPU. If True this will slow performance significantly
snapsgpu2cpu = False
# Numerical dispersion analysis parameters
# Threshold (dB) down from maximum power (0dB) of main frequency used
# to calculate highest frequency for numerical dispersion analysis
# Maximum allowable percentage physical phase-velocity phase error
# Minimum grid sampling of smallest wavelength for physical wave propagation
numdispersion = {'highestfreqthres': 40, 'maxnumericaldisp': 2, 'mingridsampling': 3}
# Simulation mode - 2D TMx, 2D TMy, 2D TMz, or 3D
mode = None
# Data types
# Solid and ID arrays use 32-bit integers (0 to 4294967295)
# Rigid arrays use 8-bit integers (the smallest available type to store true/false)
# Fractal and dispersive coefficient arrays use complex numbers (complextype)
# which are represented as two floats
# Main field arrays use floats (floattype) and complex numbers (complextype)
# Precision of floattype and complextype: single or double
precision = 'single'
if precision == 'single':
# For numpy arrays
floattype = np.float32
complextype = np.complex64
# For C (CUDA) arrays
cudafloattype = 'float'
cudacomplextype = 'pycuda::complex<float>'
elif precision == 'double':
# For numpy arrays
floattype = np.float64
complextype = np.complex128
# For C (CUDA) arrays
cudafloattype = 'double'
cudacomplextype = 'pycuda::complex<double>'

50
gprMax/constants.py
查看文件

@@ -1,50 +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 <http://www.gnu.org/licenses/>.
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
# Impedance of free space (Ohms)
z0 = np.sqrt(m0 / e0)
# Data types:
# Solid and ID arrays use 32-bit integers (0 to 4294967295)
# Rigid arrays use 8-bit integers (the smallest available type to store true/false)
# Fractal and dispersive coefficient arrays use complex numbers (complextype)
# which are represented as two floats
# Main field arrays use floats (floattype) and complex numbers (complextype)
# Single precision
# For numpy arrays
floattype = np.float32
complextype = np.complex64
# For C (CUDA) arrays
cudafloattype = 'float'
cudacomplextype = 'pycuda::complex<float>'
# Double precision
# For numpy arrays
# floattype = np.float64
# complextype = np.complex128
# For C (CUDA) arrays
# cudafloattype = 'double'
# cudacomplextype = 'pycuda::complex<double>'

4
gprMax/fields_updates_ext.pyx
查看文件

@@ -20,8 +20,8 @@ import numpy as np
cimport numpy as np
from cython.parallel import prange
from gprMax.constants cimport floattype_t
from gprMax.constants cimport complextype_t
from gprMax.config cimport floattype_t
from gprMax.config cimport complextype_t
###############################################

9
gprMax/fractals.py
查看文件

@@ -19,8 +19,9 @@
import numpy as np
from scipy import fftpack
from gprMax.constants import floattype
from gprMax.constants import complextype
from gprMax.config import floattype
from gprMax.config import complextype
from gprMax.config import hostinfo
from gprMax.fractals_generate_ext import generate_fractal2D
from gprMax.fractals_generate_ext import generate_fractal3D
from gprMax.utilities import round_value
@@ -90,7 +91,7 @@ class FractalSurface(object):
A = fftpack.fftshift(A)
# Generate fractal
generate_fractal2D(surfacedims[0], surfacedims[1], G.nthreads, self.b, self.weighting, v1, A, self.fractalsurface)
generate_fractal2D(surfacedims[0], surfacedims[1], hostinfo['ompthreads'], self.b, self.weighting, v1, A, self.fractalsurface)
# Shift the zero frequency component to start of the array
self.fractalsurface = fftpack.ifftshift(self.fractalsurface)
@@ -172,7 +173,7 @@ class FractalVolume(object):
A = fftpack.fftshift(A)
# Generate fractal
generate_fractal3D(self.nx, self.ny, self.nz, G.nthreads, self.b, self.weighting, v1, A, self.fractalvolume)
generate_fractal3D(self.nx, self.ny, self.nz, hostinfo['ompthreads'], self.b, self.weighting, v1, A, self.fractalvolume)
# Shift the zero frequency component to the start of the array
self.fractalvolume = fftpack.ifftshift(self.fractalvolume)

2
gprMax/fractals_generate_ext.pyx
查看文件

@@ -20,7 +20,7 @@ import numpy as np
cimport numpy as np
from cython.parallel import prange
from gprMax.constants cimport complextype_t
from gprMax.config cimport complextype_t
cpdef void generate_fractal2D(int nx, int ny, int nthreads, int b, np.float64_t[:] weighting, np.float64_t[:] v1, np.complex128_t[:, ::1] A, complextype_t[:, ::1] fractalsurface):

126
gprMax/gprMax.py
查看文件

@@ -28,15 +28,17 @@ from enum import Enum
import h5py
import numpy as np
from gprMax._version import __version__, codename
from gprMax.constants import c
from gprMax.constants import e0
from gprMax.constants import m0
from gprMax.constants import z0
from gprMax._version import __version__
from gprMax._version import codename
import gprMax.config
from gprMax.config import c
from gprMax.config import e0
from gprMax.config import m0
from gprMax.config import z0
from gprMax.config import hostinfo
from gprMax.exceptions import GeneralError
from gprMax.model_build_run import run_model
from gprMax.utilities import detect_check_gpus
from gprMax.utilities import get_host_info
from gprMax.utilities import get_terminal_width
from gprMax.utilities import human_size
from gprMax.utilities import logo
@@ -60,7 +62,6 @@ def main():
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')
parser.add_argument('--opt-taguchi', action='store_true', default=False, help='flag to optimise parameters using the Taguchi optimisation method')
args = parser.parse_args()
run_main(args)
@@ -79,7 +80,6 @@ def api(
geometry_only=False,
geometry_fixed=False,
write_processed=False,
opt_taguchi=False
):
"""If installed as a module this is the entry point."""
@@ -100,7 +100,6 @@ def api(
args.geometry_only = geometry_only
args.geometry_fixed = geometry_fixed
args.write_processed = write_processed
args.opt_taguchi = opt_taguchi
run_main(args)
@@ -116,30 +115,23 @@ def run_main(args):
# Print gprMax logo, version, and licencing/copyright information
logo(__version__ + ' (' + codename + ')')
# Print information about host machine
hyperthreading = ', {} cores with HT'.format(hostinfo['logicalcores']) if hostinfo['hyperthreading'] else ''
print('\nHost: {} | {} | {} x {} ({} cores{}) | {} RAM | {}'.format(hostinfo['hostname'],
hostinfo['machineID'], hostinfo['sockets'], hostinfo['cpuID'], hostinfo['physicalcores'],
hyperthreading, human_size(hostinfo['ram'], a_kilobyte_is_1024_bytes=True), hostinfo['osversion']))
# Get information/setup any Nvidia GPU(s)
if args.gpu is not None:
# Flatten a list of lists
if any(isinstance(element, list) for element in args.gpu):
args.gpu = [val for sublist in args.gpu for val in sublist]
gpus1, allgpustext = detect_check_gpus(args.gpu)
print('with GPU(s): {}'.format(' | '.join(allgpustext)))
# Process input file
with open_path_file(args.inputfile) as inputfile:
# Get information about host machine
hostinfo = get_host_info()
hyperthreading = ', {} cores with Hyper-Threading'.format(hostinfo['logicalcores']) if hostinfo['hyperthreading'] else ''
print('\nHost: {} | {} | {} x {} ({} cores{}) | {} RAM | {}'.format(hostinfo['hostname'],
hostinfo['machineID'], hostinfo['sockets'], hostinfo['cpuID'], hostinfo['physicalcores'],
hyperthreading, human_size(hostinfo['ram'], a_kilobyte_is_1024_bytes=True), hostinfo['osversion']))
# Get information/setup any Nvidia GPU(s)
if args.gpu is not None:
# Flatten a list of lists
if any(isinstance(element, list) for element in args.gpu):
args.gpu = [val for sublist in args.gpu for val in sublist]
gpus, allgpustext = detect_check_gpus(args.gpu)
print('GPU(s) detected: {}'.format(' | '.join(allgpustext)))
# If in MPI mode or benchmarking provide list of GPU objects, otherwise
# provide single GPU object
if args.mpi or args.mpi_no_spawn or args.benchmark:
args.gpu = gpus
else:
args.gpu = gpus[0]
# Create a separate namespace that users can access in any Python code blocks in the input file
usernamespace = {'c': c, 'e0': e0, 'm0': m0, 'z0': z0, 'number_model_runs': args.n, 'inputfile': os.path.abspath(inputfile.name)}
@@ -147,20 +139,10 @@ def run_main(args):
# Process for benchmarking simulation #
#######################################
if args.benchmark:
if args.mpi or args.opt_taguchi or args.task or args.n > 1:
raise GeneralError('Benchmarking mode cannot be combined with MPI, job array, or Taguchi optimisation modes, or multiple model runs.')
if args.mpi or args.task or args.n > 1:
raise GeneralError('Benchmarking mode cannot be combined with MPI, job array mode, or multiple model runs.')
run_benchmark_sim(args, inputfile, usernamespace)
####################################################
# Process for simulation with Taguchi optimisation #
####################################################
elif args.opt_taguchi:
if args.mpi_worker: # Special case for MPI spawned workers - they do not need to enter the Taguchi optimisation mode
run_mpi_sim(args, inputfile, usernamespace)
else:
from gprMax.optimisation_taguchi import run_opt_sim
run_opt_sim(args, inputfile, usernamespace)
################################################
# Process for standard simulation (CPU or GPU) #
################################################
@@ -185,10 +167,12 @@ def run_main(args):
else:
if args.task and args.restart:
raise GeneralError('Job array and restart modes cannot be used together')
if gprMax.config.gpus:
gprMax.config.gpus = gprMax.config.gpus[0]
run_std_sim(args, inputfile, usernamespace)
def run_std_sim(args, inputfile, usernamespace, optparams=None):
def run_std_sim(args, inputfile, usernamespace):
"""
Run standard simulation - models are run one after another and each model
is parallelised using either OpenMP (CPU) or CUDA (GPU)
@@ -198,8 +182,6 @@ def run_std_sim(args, inputfile, usernamespace, optparams=None):
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.
optparams (dict): Optional argument. For Taguchi optimisation it
provides the parameters to optimise and their values.
"""
# Set range for number of models to run
@@ -217,18 +199,9 @@ def run_std_sim(args, inputfile, usernamespace, optparams=None):
tsimstart = timer()
for currentmodelrun in range(modelstart, modelend):
# If Taguchi optimistaion, add specific value for each parameter to
# optimise for each experiment to user accessible namespace
if optparams:
tmp = {}
tmp.update((key, value[currentmodelrun - 1]) for key, value in optparams.items())
modelusernamespace = usernamespace.copy()
modelusernamespace.update({'optparams': tmp})
else:
modelusernamespace = usernamespace
run_model(args, currentmodelrun, modelend - 1, numbermodelruns, inputfile, modelusernamespace)
run_model(args, currentmodelrun, modelend - 1, numbermodelruns, inputfile, usernamespace)
tsimend = timer()
simcompletestr = '\n=== Simulation completed in [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=tsimend - tsimstart))
simcompletestr = '\n=== Simulation on {} completed in [HH:MM:SS]: {}'.format(hostinfo['hostname'], datetime.timedelta(seconds=tsimend - tsimstart))
print('{} {}\n'.format(simcompletestr, '=' * (get_terminal_width() - 1 - len(simcompletestr))))
@@ -245,9 +218,8 @@ def run_benchmark_sim(args, inputfile, usernamespace):
Python code blocks in input file.
"""
# Get information about host machine
hostinfo = get_host_info()
hyperthreading = ', {} cores with Hyper-Threading'.format(hostinfo['logicalcores']) if hostinfo['hyperthreading'] else ''
# Store information about host machine
hyperthreading = ', {} cores with HT'.format(hostinfo['logicalcores']) if hostinfo['hyperthreading'] else ''
machineIDlong = '{}; {} x {} ({} cores{}); {} RAM; {}'.format(hostinfo['machineID'], hostinfo['sockets'], hostinfo['cpuID'], hostinfo['physicalcores'], hyperthreading, human_size(hostinfo['ram'], a_kilobyte_is_1024_bytes=True), hostinfo['osversion'])
# Initialise arrays to hold CPU thread info and times, and GPU info and times
@@ -319,7 +291,7 @@ def run_benchmark_sim(args, inputfile, usernamespace):
print('{} {}\n'.format(simcompletestr, '=' * (get_terminal_width() - 1 - len(simcompletestr))))
def run_mpi_sim(args, inputfile, usernamespace, optparams=None):
def run_mpi_sim(args, inputfile, usernamespace):
"""
Run mixed mode MPI/OpenMP simulation - MPI task farm for models with
each model parallelised using either OpenMP (CPU) or CUDA (GPU)
@@ -329,8 +301,6 @@ def run_mpi_sim(args, inputfile, usernamespace, optparams=None):
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.
optparams (dict): Optional argument. For Taguchi optimisation it
provides the parameters to optimise and their values.
"""
from mpi4py import MPI
@@ -436,26 +406,16 @@ def run_mpi_sim(args, inputfile, usernamespace, optparams=None):
for work in iter(lambda: comm.sendrecv(0, dest=0), StopIteration):
currentmodelrun = work['currentmodelrun']
# If Taguchi optimisation, add specific value for each parameter to
# optimise for each experiment to user accessible namespace
if 'optparams' in work:
tmp = {}
tmp.update((key, value[currentmodelrun - 1]) for key, value in work['optparams'].items())
modelusernamespace = usernamespace.copy()
modelusernamespace.update({'optparams': tmp})
else:
modelusernamespace = usernamespace
# Run the model
print('Starting MPI spawned worker (parent: {}, rank: {}) on {} with model {}/{}{}\n'.format(work['mpicommname'], rank, hostname, currentmodelrun, numbermodelruns, gpuinfo))
tsolve = run_model(args, currentmodelrun, modelend - 1, numbermodelruns, inputfile, modelusernamespace)
tsolve = run_model(args, currentmodelrun, modelend - 1, numbermodelruns, inputfile, usernamespace)
print('Completed MPI spawned worker (parent: {}, rank: {}) on {} with model {}/{}{} in [HH:MM:SS]: {}\n'.format(work['mpicommname'], rank, hostname, currentmodelrun, numbermodelruns, gpuinfo, datetime.timedelta(seconds=tsolve)))
# Shutdown
comm.Disconnect()
def run_mpi_no_spawn_sim(args, inputfile, usernamespace, optparams=None):
def run_mpi_no_spawn_sim(args, inputfile, usernamespace):
"""
Alternate MPI implementation that avoids using the MPI spawn mechanism.
This implementation is designed to be used as
@@ -469,8 +429,6 @@ def run_mpi_no_spawn_sim(args, inputfile, usernamespace, optparams=None):
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.
optparams (dict): Optional argument. For Taguchi optimisation it
provides the parameters to optimise and their values.
"""
from mpi4py import MPI
@@ -547,20 +505,8 @@ def run_mpi_no_spawn_sim(args, inputfile, usernamespace, optparams=None):
# Run a model
if tag == tags.START.value:
# If Taguchi optimistaion, add specific value for each parameter
# to optimise for each experiment to user accessible namespace
if optparams:
tmp = {}
tmp.update((key, value[currentmodelrun - 1]) for key, value in optparams.items())
modelusernamespace = usernamespace.copy()
modelusernamespace.update({'optparams': tmp})
else:
modelusernamespace = usernamespace
# Run the model
print('Starting MPI worker (parent: {}, rank: {}) on {} with model {}/{}{}\n'.format(comm.name, rank, hostname, currentmodelrun, numbermodelruns, gpuinfo))
tsolve = run_model(args, currentmodelrun, modelend - 1, numbermodelruns, inputfile, modelusernamespace)
tsolve = run_model(args, currentmodelrun, modelend - 1, numbermodelruns, inputfile, usernamespace)
comm.send(None, dest=0, tag=tags.DONE.value)
print('Completed MPI worker (parent: {}, rank: {}) on {} with model {}/{}{} in [HH:MM:SS]: {}\n'.format(comm.name, rank, hostname, currentmodelrun, numbermodelruns, gpuinfo, datetime.timedelta(seconds=tsolve)))

63
gprMax/grid.py
查看文件

@@ -25,9 +25,13 @@ init()
import numpy as np
np.seterr(invalid='raise')
from gprMax.constants import c
from gprMax.constants import floattype
from gprMax.constants import complextype
import gprMax.config
from gprMax.config import c
from gprMax.config import floattype
from gprMax.config import complextype
from gprMax.config import numdispersion
from gprMax.config import gpus as gpu
from gprMax.config import hostinfo
from gprMax.exceptions import GeneralError
from gprMax.materials import Material
from gprMax.pml import PML
@@ -88,36 +92,8 @@ class FDTDGrid(Grid):
self.inputdirectory = ''
self.outputdirectory = ''
self.title = ''
self.messages = True
self.progressbars = self.messages
self.memoryusage = 0
# Get information about host machine
self.hostinfo = None
# CPU - OpenMP threads
self.nthreads = 0
# GPU
# Threads per block - electric and magnetic field updates
self.tpb = (256, 1, 1)
# GPU object
self.gpu = None
# Copy snapshot data from GPU to CPU during simulation
# N.B. This will happen if the requested snapshots are too large to fit
# on the memory of the GPU. If True this will slow performance significantly
self.snapsgpu2cpu = False
# Threshold (dB) down from maximum power (0dB) of main frequency used
# to calculate highest frequency for numerical dispersion analysis
self.highestfreqthres = 40
# Maximum allowable percentage physical phase-velocity phase error
self.maxnumericaldisp = 2
# Minimum grid sampling of smallest wavelength for physical wave propagation
self.mingridsampling = 3
self.nx = 0
self.ny = 0
self.nz = 0
@@ -125,7 +101,6 @@ class FDTDGrid(Grid):
self.dy = 0
self.dz = 0
self.dt = 0
self.mode = None
self.iterations = 0
self.timewindow = 0
@@ -232,21 +207,21 @@ class FDTDGrid(Grid):
"""
# Check if model can be built and/or run on host
if self.memoryusage > self.hostinfo['ram']:
raise GeneralError('Memory (RAM) required ~{} exceeds {} detected!\n'.format(human_size(self.memoryusage), human_size(self.hostinfo['ram'], a_kilobyte_is_1024_bytes=True)))
if self.memoryusage > hostinfo['ram']:
raise GeneralError('Memory (RAM) required ~{} exceeds {} detected!\n'.format(human_size(self.memoryusage), human_size(hostinfo['ram'], a_kilobyte_is_1024_bytes=True)))
# Check if model can be run on specified GPU if required
if self.gpu is not None:
if self.memoryusage - snapsmemsize > self.gpu.totalmem:
raise GeneralError('Memory (RAM) required ~{} exceeds {} detected on specified {} - {} GPU!\n'.format(human_size(self.memoryusage), human_size(self.gpu.totalmem, a_kilobyte_is_1024_bytes=True), self.gpu.deviceID, self.gpu.name))
if gpu is not None:
if self.memoryusage - snapsmemsize > gpu.totalmem:
raise GeneralError('Memory (RAM) required ~{} exceeds {} detected on specified {} - {} GPU!\n'.format(human_size(self.memoryusage), human_size(gpu.totalmem, a_kilobyte_is_1024_bytes=True), gpu.deviceID, gpu.name))
# If the required memory without the snapshots will fit on the GPU then transfer and store snaphots on host
if snapsmemsize != 0 and self.memoryusage - snapsmemsize < self.gpu.totalmem:
self.snapsgpu2cpu = True
if snapsmemsize != 0 and self.memoryusage - snapsmemsize < gpu.totalmem:
gprMax.config.snapsgpu2cpu = True
def gpu_set_blocks_per_grid(self):
"""Set the blocks per grid size used for updating the electric and magnetic field arrays on a GPU."""
self.bpg = (int(np.ceil(((self.nx + 1) * (self.ny + 1) * (self.nz + 1)) / self.tpb[0])), 1, 1)
gprMax.config.gpus.bpg = (int(np.ceil(((self.nx + 1) * (self.ny + 1) * (self.nz + 1)) / gpu.tpb[0])), 1, 1)
def gpu_initialise_arrays(self):
"""Initialise standard field arrays on GPU."""
@@ -324,7 +299,7 @@ def dispersion_analysis(G):
# Set maximum frequency to a threshold drop from maximum power, ignoring DC value
try:
freqthres = np.where(power[freqmaxpower:] < -G.highestfreqthres)[0][0] + freqmaxpower
freqthres = np.where(power[freqmaxpower:] < -numdispersion['highestfreqthres'])[0][0] + freqmaxpower
results['maxfreq'].append(freqs[freqthres])
except ValueError:
results['error'] = 'unable to calculate maximum power from waveform, most likely due to undersampling.'
@@ -365,9 +340,9 @@ def dispersion_analysis(G):
minwavelength = minvelocity / results['maxfreq']
# Maximum spatial step
if '3D' in G.mode:
if '3D' in gprMax.config.mode:
delta = max(G.dx, G.dy, G.dz)
elif '2D' in G.mode:
elif '2D' in gprMax.config.mode:
if G.nx == 1:
delta = max(G.dy, G.dz)
elif G.ny == 1:
@@ -382,7 +357,7 @@ def dispersion_analysis(G):
results['N'] = minwavelength / delta
# Check grid sampling will result in physical wave propagation
if int(np.floor(results['N'])) >= G.mingridsampling:
if int(np.floor(results['N'])) >= numdispersion['mingridsampling']:
# Numerical phase velocity
vp = np.pi / (results['N'] * np.arcsin((1 / S) * np.sin((np.pi * S) / results['N'])))

36
gprMax/input_cmds_geometry.py
查看文件

@@ -23,7 +23,9 @@ import h5py
import numpy as np
from tqdm import tqdm
from gprMax.constants import floattype
from gprMax.config import floattype
from gprMax.config import messages
from gprMax.config import progressbars
from gprMax.input_cmds_file import check_cmd_names
from gprMax.input_cmds_multiuse import process_multicmds
from gprMax.exceptions import CmdInputError
@@ -61,11 +63,11 @@ def process_geometrycmds(geometry, G):
# Disable progress bar if on Windows as it does not update properly
# when messages are printed for geometry
if sys.platform == 'win32':
progressbars = False
progressbarsgeo = False
else:
progressbars = not G.progressbars
progressbarsgeo = not progressbars
for object in tqdm(geometry, desc='Processing geometry related cmds', unit='cmds', ncols=get_terminal_width() - 1, file=sys.stdout, disable=progressbars):
for object in tqdm(geometry, desc='Processing geometry related cmds', unit='cmds', ncols=get_terminal_width() - 1, file=sys.stdout, disable=progressbarsgeo):
tmp = object.split()
if tmp[0] == '#geometry_objects_read:':
@@ -135,12 +137,12 @@ def process_geometrycmds(geometry, G):
G.rigidE[:, xs:xs + rigidE.shape[1], ys:ys + rigidE.shape[2], zs:zs + rigidE.shape[3]] = rigidE
G.rigidH[:, xs:xs + rigidH.shape[1], ys:ys + rigidH.shape[2], zs:zs + rigidH.shape[3]] = rigidH
G.ID[:, xs:xs + ID.shape[1], ys:ys + ID.shape[2], zs:zs + ID.shape[3]] = ID + numexistmaterials
if G.messages:
if messages:
tqdm.write('Geometry objects from file {} inserted at {:g}m, {:g}m, {:g}m, with corresponding materials file {}.'.format(geofile, xs * G.dx, ys * G.dy, zs * G.dz, matfile))
except KeyError:
averaging = False
build_voxels_from_array(xs, ys, zs, numexistmaterials, averaging, data, G.solid, G.rigidE, G.rigidH, G.ID)
if G.messages:
if messages:
tqdm.write('Geometry objects from file (voxels only) {} inserted at {:g}m, {:g}m, {:g}m, with corresponding materials file {}.'.format(geofile, xs * G.dx, ys * G.dy, zs * G.dz, matfile))
elif tmp[0] == '#edge:':
@@ -199,7 +201,7 @@ def process_geometrycmds(geometry, G):
for k in range(zs, zf):
build_edge_z(xs, ys, k, material.numID, G.rigidE, G.rigidH, G.ID)
if G.messages:
if messages:
tqdm.write('Edge from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material {} created.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, tmp[7]))
elif tmp[0] == '#plate:':
@@ -307,7 +309,7 @@ def process_geometrycmds(geometry, G):
for j in range(ys, yf):
build_face_xy(i, j, zs, numIDx, numIDy, G.rigidE, G.rigidH, G.ID)
if G.messages:
if messages:
tqdm.write('Plate from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m of material(s) {} created.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, ', '.join(materialsrequested)))
elif tmp[0] == '#triangle:':
@@ -420,7 +422,7 @@ def process_geometrycmds(geometry, G):
build_triangle(x1, y1, z1, x2, y2, z2, x3, y3, z3, normal, thickness, G.dx, G.dy, G.dz, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID)
if G.messages:
if messages:
if thickness > 0:
if averaging:
dielectricsmoothing = 'on'
@@ -515,7 +517,7 @@ def process_geometrycmds(geometry, G):
build_box(xs, xf, ys, yf, zs, zf, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID)
if G.messages:
if messages:
if averaging:
dielectricsmoothing = 'on'
else:
@@ -596,7 +598,7 @@ def process_geometrycmds(geometry, G):
build_cylinder(x1, y1, z1, x2, y2, z2, r, G.dx, G.dy, G.dz, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID)
if G.messages:
if messages:
if averaging:
dielectricsmoothing = 'on'
else:
@@ -717,7 +719,7 @@ def process_geometrycmds(geometry, G):
build_cylindrical_sector(ctr1, ctr2, level, sectorstartangle, sectorangle, r, normal, thickness, G.dx, G.dy, G.dz, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID)
if G.messages:
if messages:
if thickness > 0:
if averaging:
dielectricsmoothing = 'on'
@@ -796,7 +798,7 @@ def process_geometrycmds(geometry, G):
build_sphere(xc, yc, zc, r, G.dx, G.dy, G.dz, numID, numIDx, numIDy, numIDz, averaging, G.solid, G.rigidE, G.rigidH, G.ID)
if G.messages:
if messages:
if averaging:
dielectricsmoothing = 'on'
else:
@@ -877,7 +879,7 @@ def process_geometrycmds(geometry, G):
volume.weighting = np.array([float(tmp[8]), float(tmp[9]), float(tmp[10])])
volume.averaging = averagefractalbox
if G.messages:
if messages:
if volume.averaging:
dielectricsmoothing = 'on'
else:
@@ -1000,7 +1002,7 @@ def process_geometrycmds(geometry, G):
surface.generate_fractal_surface(G)
volume.fractalsurfaces.append(surface)
if G.messages:
if messages:
tqdm.write('Fractal surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with fractal dimension {:g}, fractal weightings {:g}, {:g}, fractal seeding {}, and range {:g}m to {:g}m, added to {}.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, surface.dimension, surface.weighting[0], surface.weighting[1], surface.seed, float(tmp[10]), float(tmp[11]), surface.operatingonID))
if tmp[0] == '#add_surface_water:':
@@ -1108,7 +1110,7 @@ def process_geometrycmds(geometry, G):
if testwater:
raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires the time step for the model to be less than the relaxation time required to model water.')
if G.messages:
if messages:
tqdm.write('Water on surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with depth {:g}m, added to {}.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, filldepth, surface.operatingonID))
if tmp[0] == '#add_grass:':
@@ -1258,7 +1260,7 @@ def process_geometrycmds(geometry, G):
volume.fractalsurfaces.append(surface)
if G.messages:
if messages:
tqdm.write('{} blades of grass on surface from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with fractal dimension {:g}, fractal seeding {}, and range {:g}m to {:g}m, added to {}.'.format(numblades, xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, surface.dimension, surface.seed, float(tmp[8]), float(tmp[9]), surface.operatingonID))
# Process any modifications to the original fractal box then generate it

72
gprMax/input_cmds_multiuse.py
查看文件

@@ -23,8 +23,12 @@ init()
import numpy as np
from tqdm import tqdm
from gprMax.constants import z0
from gprMax.constants import floattype
import gprMax.config
from gprMax.config import z0
from gprMax.config import floattype
from gprMax.config import gpus
from gprMax.config import messages
from gprMax.exceptions import CmdInputError
from gprMax.geometry_outputs import GeometryView
from gprMax.geometry_outputs import GeometryObjects
@@ -80,7 +84,7 @@ def process_multicmds(multicmds, G):
w.amp = float(tmp[1])
w.freq = float(tmp[2])
if G.messages:
if messages:
print('Waveform {} of type {} with maximum amplitude scaling {:g}, frequency {:g}Hz created.'.format(w.ID, w.type, w.amp, w.freq))
G.waveforms.append(w)
@@ -97,11 +101,11 @@ def process_multicmds(multicmds, G):
polarisation = tmp[0].lower()
if polarisation not in ('x', 'y', 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
if '2D TMx' in G.mode and (polarisation == 'y' or polarisation == 'z'):
if '2D TMx' in gprMax.config.mode and (polarisation == 'y' or polarisation == 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x in 2D TMx mode')
elif '2D TMy' in G.mode and (polarisation == 'x' or polarisation == 'z'):
elif '2D TMy' in gprMax.config.mode and (polarisation == 'x' or polarisation == 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be y in 2D TMy mode')
elif '2D TMz' in G.mode and (polarisation == 'x' or polarisation == 'y'):
elif '2D TMz' in gprMax.config.mode and (polarisation == 'x' or polarisation == 'y'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be z in 2D TMz mode')
xcoord = G.calculate_coord('x', tmp[1])
@@ -151,7 +155,7 @@ def process_multicmds(multicmds, G):
v.calculate_waveform_values(G)
if G.messages:
if messages:
print('Voltage source with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(v.polarisation, v.xcoord * G.dx, v.ycoord * G.dy, v.zcoord * G.dz, v.resistance) + startstop + 'using waveform {} created.'.format(v.waveformID))
G.voltagesources.append(v)
@@ -168,11 +172,11 @@ def process_multicmds(multicmds, G):
polarisation = tmp[0].lower()
if polarisation not in ('x', 'y', 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
if '2D TMx' in G.mode and (polarisation == 'y' or polarisation == 'z'):
if '2D TMx' in gprMax.config.mode and (polarisation == 'y' or polarisation == 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x in 2D TMx mode')
elif '2D TMy' in G.mode and (polarisation == 'x' or polarisation == 'z'):
elif '2D TMy' in gprMax.config.mode and (polarisation == 'x' or polarisation == 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be y in 2D TMy mode')
elif '2D TMz' in G.mode and (polarisation == 'x' or polarisation == 'y'):
elif '2D TMz' in gprMax.config.mode and (polarisation == 'x' or polarisation == 'y'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be z in 2D TMz mode')
xcoord = G.calculate_coord('x', tmp[1])
@@ -229,8 +233,8 @@ def process_multicmds(multicmds, G):
h.calculate_waveform_values(G)
if G.messages:
if G.mode == '2D':
if messages:
if '2D' in gprMax.config.mode:
print('Hertzian dipole is a line source in 2D with polarity {} at {:g}m, {:g}m, {:g}m,'.format(h.polarisation, h.xcoord * G.dx, h.ycoord * G.dy, h.zcoord * G.dz) + startstop + 'using waveform {} created.'.format(h.waveformID))
else:
print('Hertzian dipole with polarity {} at {:g}m, {:g}m, {:g}m,'.format(h.polarisation, h.xcoord * G.dx, h.ycoord * G.dy, h.zcoord * G.dz) + startstop + 'using waveform {} created.'.format(h.waveformID))
@@ -249,11 +253,11 @@ def process_multicmds(multicmds, G):
polarisation = tmp[0].lower()
if polarisation not in ('x', 'y', 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
if '2D TMx' in G.mode and (polarisation == 'y' or polarisation == 'z'):
if '2D TMx' in gprMax.config.mode and (polarisation == 'y' or polarisation == 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x in 2D TMx mode')
elif '2D TMy' in G.mode and (polarisation == 'x' or polarisation == 'z'):
elif '2D TMy' in gprMax.config.mode and (polarisation == 'x' or polarisation == 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be y in 2D TMy mode')
elif '2D TMz' in G.mode and (polarisation == 'x' or polarisation == 'y'):
elif '2D TMz' in gprMax.config.mode and (polarisation == 'x' or polarisation == 'y'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be z in 2D TMz mode')
xcoord = G.calculate_coord('x', tmp[1])
@@ -301,7 +305,7 @@ def process_multicmds(multicmds, G):
m.calculate_waveform_values(G)
if G.messages:
if messages:
print('Magnetic dipole with polarity {} at {:g}m, {:g}m, {:g}m,'.format(m.polarisation, m.xcoord * G.dx, m.ycoord * G.dy, m.zcoord * G.dz) + startstop + 'using waveform {} created.'.format(m.waveformID))
G.magneticdipoles.append(m)
@@ -315,18 +319,18 @@ def process_multicmds(multicmds, G):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' requires at least six parameters')
# Warn about using a transmission line on GPU
if G.gpu is not None:
if gpus is not None:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' A #transmission_line cannot currently be used with GPU solving. Consider using a #voltage_source instead.')
# Check polarity & position parameters
polarisation = tmp[0].lower()
if polarisation not in ('x', 'y', 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x, y, or z')
if '2D TMx' in G.mode and (polarisation == 'y' or polarisation == 'z'):
if '2D TMx' in gprMax.config.mode and (polarisation == 'y' or polarisation == 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be x in 2D TMx mode')
elif '2D TMy' in G.mode and (polarisation == 'x' or polarisation == 'z'):
elif '2D TMy' in gprMax.config.mode and (polarisation == 'x' or polarisation == 'z'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be y in 2D TMy mode')
elif '2D TMz' in G.mode and (polarisation == 'x' or polarisation == 'y'):
elif '2D TMz' in gprMax.config.mode and (polarisation == 'x' or polarisation == 'y'):
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' polarisation must be z in 2D TMz mode')
xcoord = G.calculate_coord('x', tmp[1])
@@ -377,7 +381,7 @@ def process_multicmds(multicmds, G):
t.calculate_waveform_values(G)
t.calculate_incident_V_I(G)
if G.messages:
if messages:
print('Transmission line with polarity {} at {:g}m, {:g}m, {:g}m, resistance {:.1f} Ohms,'.format(t.polarisation, t.xcoord * G.dx, t.ycoord * G.dy, t.zcoord * G.dz, t.resistance) + startstop + 'using waveform {} created.'.format(t.waveformID))
G.transmissionlines.append(t)
@@ -414,7 +418,7 @@ def process_multicmds(multicmds, G):
else:
r.ID = tmp[3]
# Get allowable outputs
if G.gpu is not None:
if gpus is not None:
allowableoutputs = Rx.gpu_allowableoutputs
else:
allowableoutputs = Rx.allowableoutputs
@@ -425,7 +429,7 @@ def process_multicmds(multicmds, G):
else:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' contains an output type that is not allowable. Allowable outputs in current context are {}'.format(allowableoutputs))
if G.messages:
if messages:
print('Receiver at {:g}m, {:g}m, {:g}m with output component(s) {} created.'.format(r.xcoord * G.dx, r.ycoord * G.dy, r.zcoord * G.dz, ', '.join(r.outputs)))
G.rxs.append(r)
@@ -475,7 +479,7 @@ def process_multicmds(multicmds, G):
else:
raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the step size should not be less than the spatial discretisation')
if G.messages:
if messages:
print('Receiver array {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m with steps {:g}m, {:g}m, {:g}m'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, dx * G.dx, dy * G.dy, dz * G.dz))
for x in range(xs, xf + 1, dx):
@@ -491,7 +495,7 @@ def process_multicmds(multicmds, G):
r.ID = r.__class__.__name__ + '(' + str(x) + ',' + str(y) + ',' + str(z) + ')'
for key in Rx.defaultoutputs:
r.outputs[key] = np.zeros(G.iterations, dtype=floattype)
if G.messages:
if messages:
print(' Receiver at {:g}m, {:g}m, {:g}m with output component(s) {} created.'.format(r.xcoord * G.dx, r.ycoord * G.dy, r.zcoord * G.dz, ', '.join(r.outputs)))
G.rxs.append(r)
@@ -540,7 +544,7 @@ def process_multicmds(multicmds, G):
s = Snapshot(xs, ys, zs, xf, yf, zf, dx, dy, dz, time, tmp[10])
if G.messages:
if messages:
print('Snapshot from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, at {:g} secs with filename {} created.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, dx * G.dx, dy * G.dy, dz * G.dz, s.time * G.dt, s.basefilename))
G.snapshots.append(s)
@@ -578,7 +582,7 @@ def process_multicmds(multicmds, G):
if m.se == float('inf'):
m.averagable = False
if G.messages:
if messages:
tqdm.write('Material {} with eps_r={:g}, sigma={:g} S/m; mu_r={:g}, sigma*={:g} Ohm/m created.'.format(m.ID, m.er, m.se, m.mr, m.sm))
# Append the new material object to the materials list
@@ -617,7 +621,7 @@ def process_multicmds(multicmds, G):
if material.poles > Material.maxpoles:
Material.maxpoles = material.poles
if G.messages:
if messages:
tqdm.write('Debye disperion added to {} with delta_eps_r={}, and tau={} secs created.'.format(material.ID, ', '.join('%4.2f' % deltaer for deltaer in material.deltaer), ', '.join('%4.3e' % tau for tau in material.tau)))
cmdname = '#add_dispersion_lorentz'
@@ -653,7 +657,7 @@ def process_multicmds(multicmds, G):
if material.poles > Material.maxpoles:
Material.maxpoles = material.poles
if G.messages:
if messages:
tqdm.write('Lorentz disperion added to {} with delta_eps_r={}, omega={} secs, and gamma={} created.'.format(material.ID, ', '.join('%4.2f' % deltaer for deltaer in material.deltaer), ', '.join('%4.3e' % tau for tau in material.tau), ', '.join('%4.3e' % alpha for alpha in material.alpha)))
cmdname = '#add_dispersion_drude'
@@ -688,7 +692,7 @@ def process_multicmds(multicmds, G):
if material.poles > Material.maxpoles:
Material.maxpoles = material.poles
if G.messages:
if messages:
tqdm.write('Drude disperion added to {} with omega={} secs, and gamma={} secs created.'.format(material.ID, ', '.join('%4.3e' % tau for tau in material.tau), ', '.join('%4.3e' % alpha for alpha in material.alpha)))
cmdname = '#soil_peplinski'
@@ -715,7 +719,7 @@ def process_multicmds(multicmds, G):
# Create a new instance of the Material class material (start index after pec & free_space)
s = PeplinskiSoil(tmp[6], float(tmp[0]), float(tmp[1]), float(tmp[2]), float(tmp[3]), (float(tmp[4]), float(tmp[5])))
if G.messages:
if messages:
print('Mixing model (Peplinski) used to create {} with sand fraction {:g}, clay fraction {:g}, bulk density {:g}g/cm3, sand particle density {:g}g/cm3, and water volumetric fraction {:g} to {:g} created.'.format(s.ID, s.S, s.C, s.rb, s.rs, s.mu[0], s.mu[1]))
# Append the new material object to the materials list
@@ -765,7 +769,7 @@ def process_multicmds(multicmds, G):
g = GeometryView(xs, ys, zs, xf, yf, zf, dx, dy, dz, tmp[9], fileext)
if G.messages:
if messages:
print('Geometry view from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, with filename base {} created.'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, dx * G.dx, dy * G.dy, dz * G.dz, g.basefilename))
# Append the new GeometryView object to the geometry views list
@@ -795,7 +799,7 @@ def process_multicmds(multicmds, G):
g = GeometryObjects(xs, ys, zs, xf, yf, zf, tmp[6])
if G.messages:
if messages:
print('Geometry objects in the volume from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, will be written to {}, with materials written to {}'.format(xs * G.dx, ys * G.dy, zs * G.dz, xf * G.dx, yf * G.dy, zf * G.dz, g.filename, g.materialsfilename))
# Append the new GeometryView object to the geometry objects to write list
@@ -845,7 +849,7 @@ def process_multicmds(multicmds, G):
cfs.kappa = cfskappa
cfs.sigma = cfssigma
if G.messages:
if messages:
print('PML CFS parameters: alpha (scaling: {}, scaling direction: {}, min: {:g}, max: {:g}), kappa (scaling: {}, scaling direction: {}, min: {:g}, max: {:g}), sigma (scaling: {}, scaling direction: {}, min: {:g}, max: {}) created.'.format(cfsalpha.scalingprofile, cfsalpha.scalingdirection, cfsalpha.min, cfsalpha.max, cfskappa.scalingprofile, cfskappa.scalingdirection, cfskappa.min, cfskappa.max, cfssigma.scalingprofile, cfssigma.scalingdirection, cfssigma.min, cfssigma.max))
G.cfs.append(cfs)

87
gprMax/input_cmds_singleuse.py
查看文件

@@ -28,12 +28,14 @@ init()
import numpy as np
from scipy import interpolate
from gprMax.constants import c
from gprMax.constants import floattype
import gprMax.config
from gprMax.config import c
from gprMax.config import floattype
from gprMax.config import gpus as gpu
from gprMax.config import hostinfo
from gprMax.exceptions import CmdInputError
from gprMax.exceptions import GeneralError
from gprMax.pml import PML
from gprMax.utilities import get_host_info
from gprMax.utilities import human_size
from gprMax.utilities import round_value
from gprMax.waveforms import Waveform
@@ -55,22 +57,12 @@ def process_singlecmds(singlecmds, G):
if len(tmp) != 1:
raise CmdInputError(cmd + ' requires exactly one parameter')
if singlecmds[cmd].lower() == 'y':
G.messages = True
gprMax.config.messages = True
elif singlecmds[cmd].lower() == 'n':
G.messages = False
gprMax.config.messages = False
else:
raise CmdInputError(cmd + ' requires input values of either y or n')
# Title
cmd = '#title'
if singlecmds[cmd] is not None:
G.title = singlecmds[cmd]
if G.messages:
print('Model title: {}'.format(G.title))
# Get information about host machine
hostinfo = get_host_info()
# Number of threads (OpenMP) to use
cmd = '#num_threads'
if sys.platform == 'darwin':
@@ -81,7 +73,7 @@ def process_singlecmds(singlecmds, G):
# os.environ['OMP_DISPLAY_ENV'] = 'TRUE' # Prints OMP version and environment variables (useful for debug)
# Catch bug with Windows Subsystem for Linux (https://github.com/Microsoft/BashOnWindows/issues/785)
if 'Microsoft' in G.hostinfo['osversion']:
if 'Microsoft' in hostinfo['osversion']:
os.environ['KMP_AFFINITY'] = 'disabled'
del os.environ['OMP_PLACES']
del os.environ['OMP_PROC_BIND']
@@ -92,24 +84,35 @@ def process_singlecmds(singlecmds, G):
raise CmdInputError(cmd + ' requires exactly one parameter to specify the number of threads to use')
if tmp[0] < 1:
raise CmdInputError(cmd + ' requires the value to be an integer not less than one')
G.nthreads = tmp[0]
os.environ['OMP_NUM_THREADS'] = str(G.nthreads)
gprMax.config.hostinfo['ompthreads'] = tmp[0]
os.environ['OMP_NUM_THREADS'] = str(gprMax.config.hostinfo['ompthreads'])
elif os.environ.get('OMP_NUM_THREADS'):
G.nthreads = int(os.environ.get('OMP_NUM_THREADS'))
gprMax.config.hostinfo['ompthreads'] = int(os.environ.get('OMP_NUM_THREADS'))
else:
# Set number of threads to number of physical CPU cores
G.nthreads = hostinfo['physicalcores']
os.environ['OMP_NUM_THREADS'] = str(G.nthreads)
gprMax.config.hostinfo['ompthreads'] = hostinfo['physicalcores']
os.environ['OMP_NUM_THREADS'] = str(gprMax.config.hostinfo['ompthreads'])
if G.messages:
print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
if G.nthreads > G.hostinfo['physicalcores']:
print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, hostinfo['physicalcores']) + Style.RESET_ALL)
if gprMax.config.messages:
print('CPU (OpenMP) threads: {}'.format(gprMax.config.hostinfo['ompthreads']))
if gprMax.config.hostinfo['ompthreads'] > hostinfo['physicalcores']:
print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(gprMax.config.hostinfo['ompthreads'], hostinfo['physicalcores']) + Style.RESET_ALL)
# Print information about any GPU in use
if G.messages:
if G.gpu is not None:
print('GPU solving using: {} - {}'.format(G.gpu.deviceID, G.gpu.name))
if gprMax.config.messages:
if gpu is not None:
print('GPU: {} - {}'.format(gpu.deviceID, gpu.name))
# Print information about precision of main field values
if gprMax.config.messages:
print('Output data type: {}\n'.format(np.dtype(floattype).name))
# Title
cmd = '#title'
if singlecmds[cmd] is not None:
G.title = singlecmds[cmd]
if gprMax.config.messages:
print('Model title: {}'.format(G.title))
# Spatial discretisation
cmd = '#dx_dy_dz'
@@ -125,7 +128,7 @@ def process_singlecmds(singlecmds, G):
G.dx = tmp[0]
G.dy = tmp[1]
G.dz = tmp[2]
if G.messages:
if gprMax.config.messages:
print('Spatial discretisation: {:g} x {:g} x {:g}m'.format(G.dx, G.dy, G.dz))
# Domain
@@ -138,35 +141,35 @@ def process_singlecmds(singlecmds, G):
G.nz = round_value(tmp[2] / G.dz)
if G.nx == 0 or G.ny == 0 or G.nz == 0:
raise CmdInputError(cmd + ' requires at least one cell in every dimension')
if G.messages:
if gprMax.config.messages:
print('Domain size: {:g} x {:g} x {:g}m ({:d} x {:d} x {:d} = {:g} cells)'.format(tmp[0], tmp[1], tmp[2], G.nx, G.ny, G.nz, (G.nx * G.ny * G.nz)))
# Time step CFL limit (either 2D or 3D); switch off appropriate PMLs for 2D
if G.nx == 1:
G.dt = 1 / (c * np.sqrt((1 / G.dy) * (1 / G.dy) + (1 / G.dz) * (1 / G.dz)))
G.mode = '2D TMx'
gprMax.config.mode = '2D TMx'
G.pmlthickness['x0'] = 0
G.pmlthickness['xmax'] = 0
elif G.ny == 1:
G.dt = 1 / (c * np.sqrt((1 / G.dx) * (1 / G.dx) + (1 / G.dz) * (1 / G.dz)))
G.mode = '2D TMy'
gprMax.config.mode = '2D TMy'
G.pmlthickness['y0'] = 0
G.pmlthickness['ymax'] = 0
elif G.nz == 1:
G.dt = 1 / (c * np.sqrt((1 / G.dx) * (1 / G.dx) + (1 / G.dy) * (1 / G.dy)))
G.mode = '2D TMz'
gprMax.config.mode = '2D TMz'
G.pmlthickness['z0'] = 0
G.pmlthickness['zmax'] = 0
else:
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'
gprMax.config.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.
G.dt = round_value(G.dt, decimalplaces=d.getcontext().prec - 1)
if G.messages:
print('Mode: {}'.format(G.mode))
if gprMax.config.messages:
print('Mode: {}'.format(gprMax.config.mode))
print('Time step (at CFL limit): {:g} secs'.format(G.dt))
# Time step stability factor
@@ -178,7 +181,7 @@ def process_singlecmds(singlecmds, G):
if tmp[0] <= 0 or tmp[0] > 1:
raise CmdInputError(cmd + ' requires the value of the time step stability factor to be between zero and one')
G.dt = G.dt * tmp[0]
if G.messages:
if gprMax.config.messages:
print('Time step (modified): {:g} secs'.format(G.dt))
# Time window
@@ -203,7 +206,7 @@ def process_singlecmds(singlecmds, G):
G.iterations = int(np.ceil(tmp / G.dt)) + 1
else:
raise CmdInputError(cmd + ' must have a value greater than zero')
if G.messages:
if gprMax.config.messages:
print('Time window: {:g} secs ({} iterations)'.format(G.timewindow, G.iterations))
# PML cells
@@ -245,7 +248,7 @@ def process_singlecmds(singlecmds, G):
G.srcsteps[0] = round_value(float(tmp[0]) / G.dx)
G.srcsteps[1] = round_value(float(tmp[1]) / G.dy)
G.srcsteps[2] = round_value(float(tmp[2]) / G.dz)
if G.messages:
if gprMax.config.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))
# rx_steps
@@ -257,7 +260,7 @@ def process_singlecmds(singlecmds, G):
G.rxsteps[0] = round_value(float(tmp[0]) / G.dx)
G.rxsteps[1] = round_value(float(tmp[1]) / G.dy)
G.rxsteps[2] = round_value(float(tmp[2]) / G.dz)
if G.messages:
if gprMax.config.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))
# Excitation file for user-defined source waveforms
@@ -281,7 +284,7 @@ def process_singlecmds(singlecmds, G):
if not os.path.isfile(excitationfile):
excitationfile = os.path.abspath(os.path.join(G.inputdirectory, excitationfile))
if G.messages:
if gprMax.config.messages:
print('\nExcitation file: {}'.format(excitationfile))
# Get waveform names
@@ -323,7 +326,7 @@ def process_singlecmds(singlecmds, G):
# Interpolate waveform values
w.userfunc = interpolate.interp1d(waveformtime, singlewaveformvalues, **kwargs)
if G.messages:
if gprMax.config.messages:
print('User waveform {} created using {} and, if required, interpolation parameters (kind: {}, fill value: {}).'.format(w.ID, timestr, kwargs['kind'], kwargs['fill_value']))
G.waveforms.append(w)

6
gprMax/materials.py
查看文件

@@ -18,9 +18,9 @@
import numpy as np
from gprMax.constants import e0
from gprMax.constants import m0
from gprMax.constants import complextype
from gprMax.config import e0
from gprMax.config import m0
from gprMax.config import complextype
class Material(object):

140
gprMax/model_build_run.py
查看文件

@@ -28,13 +28,19 @@ from colorama import Fore
from colorama import Style
init()
import numpy as np
from terminaltables import AsciiTable
from terminaltables import SingleTable
from tqdm import tqdm
from gprMax.constants import floattype
from gprMax.constants import complextype
from gprMax.constants import cudafloattype
from gprMax.constants import cudacomplextype
from gprMax.config import floattype
from gprMax.config import complextype
from gprMax.config import cudafloattype
from gprMax.config import cudacomplextype
from gprMax.config import numdispersion
from gprMax.config import hostinfo
from gprMax.config import gpus as gpu
from gprMax.config import messages
from gprMax.config import progressbars
from gprMax.config import snapsgpu2cpu
from gprMax.exceptions import GeneralError
from gprMax.fields_outputs import store_outputs
@@ -56,8 +62,8 @@ from gprMax.input_cmds_geometry import process_geometrycmds
from gprMax.input_cmds_file import process_python_include_code
from gprMax.input_cmds_file import write_processed_file
from gprMax.input_cmds_file import check_cmd_names
from gprMax.input_cmds_multiuse import process_multicmds
from gprMax.input_cmds_singleuse import process_singlecmds
from gprMax.input_cmds_multiuse import process_multicmds
from gprMax.materials import Material
from gprMax.materials import process_materials
from gprMax.pml import CFS
@@ -71,7 +77,6 @@ from gprMax.snapshots import gpu_get_snapshot_array
from gprMax.snapshots_gpu import kernel_template_store_snapshot
from gprMax.sources import gpu_initialise_src_arrays
from gprMax.source_updates_gpu import kernels_template_sources
from gprMax.utilities import get_host_info
from gprMax.utilities import get_terminal_width
from gprMax.utilities import human_size
from gprMax.utilities import open_path_file
@@ -112,18 +117,10 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
# Initialise an instance of the FDTDGrid class
G = FDTDGrid()
# Get information about host machine
# (need to save this info to FDTDGrid instance after it has been created)
G.hostinfo = get_host_info()
# Single GPU object
if args.gpu:
G.gpu = args.gpu
G.inputfilename = os.path.split(inputfile.name)[1]
G.inputdirectory = os.path.dirname(os.path.abspath(inputfile.name))
inputfilestr = '\n--- Model {}/{}, input file: {}'.format(currentmodelrun, modelend, inputfile.name)
if G.messages:
if messages:
print(Fore.GREEN + '{} {}\n'.format(inputfilestr, '-' * (get_terminal_width() - 1 - len(inputfilestr))) + Style.RESET_ALL)
# Add the current model run to namespace that can be accessed by
@@ -138,7 +135,7 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
for key, value in sorted(usernamespace.items()):
if key != '__builtins__':
uservars += '{}: {}, '.format(key, value)
if G.messages:
if messages:
print('Constants/variables used/available for Python scripting: {{{}}}\n'.format(uservars[:-2]))
# Write a file containing the input commands after Python or include file commands have been processed
@@ -160,19 +157,19 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
# Process parameters for commands that can only occur once in the model
process_singlecmds(singlecmds, G)
from gprMax.config import mode
# Process parameters for commands that can occur multiple times in the model
if G.messages: print()
if messages: print()
process_multicmds(multicmds, G)
# Estimate and check memory (RAM) usage
G.memory_estimate_basic()
G.memory_check()
if G.messages:
if G.gpu is None:
print('\nMemory (RAM) required: ~{}\n'.format(human_size(G.memoryusage)))
else:
print('\nMemory (RAM) required: ~{} host + ~{} GPU\n'.format(human_size(G.memoryusage), human_size(G.memoryusage)))
if messages:
print('\nMemory (RAM) required: ~{}\n'.format(human_size(G.memoryusage)))
if gpu:
print('\nGPU memory (RAM) required: ~{}\n'.format(human_size(G.memoryusage)))
# Initialise an array for volumetric material IDs (solid), boolean
# arrays for specifying materials not to be averaged (rigid),
@@ -186,16 +183,16 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
process_geometrycmds(geometry, G)
# Build the PMLs and calculate initial coefficients
if G.messages: print()
if messages: print()
if all(value == 0 for value in G.pmlthickness.values()):
if G.messages:
if messages:
print('PML: switched off')
pass # If all the PMLs are switched off don't need to build anything
else:
# Set default CFS parameters for PML if not given
if not G.cfs:
G.cfs = [CFS()]
if G.messages:
if messages:
if all(value == G.pmlthickness['x0'] for value in G.pmlthickness.values()):
pmlinfo = str(G.pmlthickness['x0'])
else:
@@ -204,14 +201,14 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
pmlinfo += '{}: {}, '.format(key, value)
pmlinfo = pmlinfo[:-2] + ' cells'
print('PML: formulation: {}, order: {}, thickness: {}'.format(G.pmlformulation, len(G.cfs), pmlinfo))
pbar = tqdm(total=sum(1 for value in G.pmlthickness.values() if value > 0), desc='Building PML boundaries', ncols=get_terminal_width() - 1, file=sys.stdout, disable=not G.progressbars)
pbar = tqdm(total=sum(1 for value in G.pmlthickness.values() if value > 0), desc='Building PML boundaries', ncols=get_terminal_width() - 1, file=sys.stdout, disable=not progressbars)
build_pmls(G, pbar)
pbar.close()
# Build the model, i.e. set the material properties (ID) for every edge
# of every Yee cell
if G.messages: print()
pbar = tqdm(total=2, desc='Building main grid', ncols=get_terminal_width() - 1, file=sys.stdout, disable=not G.progressbars)
if messages: print()
pbar = tqdm(total=2, desc='Building main grid', ncols=get_terminal_width() - 1, file=sys.stdout, disable=not progressbars)
build_electric_components(G.solid, G.rigidE, G.ID, G)
pbar.update()
build_magnetic_components(G.solid, G.rigidH, G.ID, G)
@@ -220,19 +217,19 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
# Add PEC boundaries to invariant direction in 2D modes
# N.B. 2D modes are a single cell slice of 3D grid
if '2D TMx' in G.mode:
if '2D TMx' in mode:
# Ey & Ez components
G.ID[1, 0, :, :] = 0
G.ID[1, 1, :, :] = 0
G.ID[2, 0, :, :] = 0
G.ID[2, 1, :, :] = 0
elif '2D TMy' in G.mode:
elif '2D TMy' in mode:
# Ex & Ez components
G.ID[0, :, 0, :] = 0
G.ID[0, :, 1, :] = 0
G.ID[2, :, 0, :] = 0
G.ID[2, :, 1, :] = 0
elif '2D TMz' in G.mode:
elif '2D TMz' in mode:
# Ex & Ey components
G.ID[0, :, :, 0] = 0
G.ID[0, :, :, 1] = 0
@@ -253,7 +250,7 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
# Update estimated memory (RAM) usage
G.memoryusage += int(3 * Material.maxpoles * (G.nx + 1) * (G.ny + 1) * (G.nz + 1) * np.dtype(complextype).itemsize)
G.memory_check()
if G.messages:
if messages:
print('\nMemory (RAM) required - updated (dispersive): ~{}\n'.format(human_size(G.memoryusage)))
G.initialise_dispersive_arrays()
@@ -266,34 +263,34 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
snapsmemsize += (2 * snap.datasizefield)
G.memoryusage += int(snapsmemsize)
G.memory_check(snapsmemsize=int(snapsmemsize))
if G.messages:
if messages:
print('\nMemory (RAM) required - updated (snapshots): ~{}\n'.format(human_size(G.memoryusage)))
# Process complete list of materials - calculate update coefficients,
# store in arrays, and build text list of materials/properties
materialsdata = process_materials(G)
if G.messages:
if messages:
print('\nMaterials:')
materialstable = AsciiTable(materialsdata)
materialstable = SingleTable(materialsdata)
materialstable.outer_border = False
materialstable.justify_columns[0] = 'right'
print(materialstable.table)
# Check to see if numerical dispersion might be a problem
results = dispersion_analysis(G)
if results['error'] and G.messages:
if results['error'] and messages:
print(Fore.RED + "\nWARNING: Numerical dispersion analysis not carried out as {}".format(results['error']) + Style.RESET_ALL)
elif results['N'] < G.mingridsampling:
elif results['N'] < numdispersion['mingridsampling']:
raise GeneralError("Non-physical wave propagation: Material '{}' has wavelength sampled by {} cells, less than required minimum for physical wave propagation. Maximum significant frequency estimated as {:g}Hz".format(results['material'].ID, results['N'], results['maxfreq']))
elif results['deltavp'] and np.abs(results['deltavp']) > G.maxnumericaldisp and G.messages:
elif results['deltavp'] and np.abs(results['deltavp']) > numdispersion['maxnumericaldisp'] and messages:
print(Fore.RED + "\nWARNING: Potentially significant numerical dispersion. 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']) + Style.RESET_ALL)
elif results['deltavp'] and G.messages:
elif results['deltavp'] and messages:
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']))
# If geometry information to be reused between model runs
else:
inputfilestr = '\n--- Model {}/{}, input file (not re-processed, i.e. geometry fixed): {}'.format(currentmodelrun, modelend, inputfile.name)
if G.messages:
if messages:
print(Fore.GREEN + '{} {}\n'.format(inputfilestr, '-' * (get_terminal_width() - 1 - len(inputfilestr))) + Style.RESET_ALL)
# Clear arrays for field components
@@ -322,18 +319,18 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
receiver.zcoord = receiver.zcoordorigin + (currentmodelrun - 1) * G.rxsteps[2]
# Write files for any geometry views and geometry object outputs
if not (G.geometryviews or G.geometryobjectswrite) and args.geometry_only and G.messages:
if not (G.geometryviews or G.geometryobjectswrite) and args.geometry_only and messages:
print(Fore.RED + '\nWARNING: No geometry views or geometry objects to output found.' + Style.RESET_ALL)
if G.geometryviews:
if G.messages: print()
if messages: print()
for i, geometryview in enumerate(G.geometryviews):
geometryview.set_filename(appendmodelnumber, G)
pbar = tqdm(total=geometryview.datawritesize, unit='byte', unit_scale=True, desc='Writing geometry view file {}/{}, {}'.format(i + 1, len(G.geometryviews), os.path.split(geometryview.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not G.progressbars)
pbar = tqdm(total=geometryview.datawritesize, unit='byte', unit_scale=True, desc='Writing geometry view file {}/{}, {}'.format(i + 1, len(G.geometryviews), os.path.split(geometryview.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not progressbars)
geometryview.write_vtk(G, pbar)
pbar.close()
if G.geometryobjectswrite:
for i, geometryobject in enumerate(G.geometryobjectswrite):
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 G.progressbars)
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 progressbars)
geometryobject.write_hdf5(G, pbar)
pbar.close()
@@ -355,17 +352,17 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
outputdir = os.path.abspath(outputdir)
if not os.path.isdir(outputdir):
os.mkdir(outputdir)
if G.messages:
if messages:
print('\nCreated output directory: {}'.format(outputdir))
# Restore current directory
os.chdir(curdir)
basename, ext = os.path.splitext(inputfilename)
outputfile = os.path.join(outputdir, basename + appendmodelnumber + '.out')
if G.messages:
if messages:
print('\nOutput file: {}\n'.format(outputfile))
# Main FDTD solving functions for either CPU or GPU
if G.gpu is None:
if gpu is None:
tsolve = solve_cpu(currentmodelrun, modelend, G)
else:
tsolve, memsolve = solve_gpu(currentmodelrun, modelend, G)
@@ -380,19 +377,18 @@ def run_model(args, currentmodelrun, modelend, numbermodelruns, inputfile, usern
if not os.path.exists(snapshotdir):
os.mkdir(snapshotdir)
if G.messages: print()
if messages: print()
for i, snap in enumerate(G.snapshots):
snap.filename = os.path.abspath(os.path.join(snapshotdir, snap.basefilename + '.vti'))
pbar = tqdm(total=snap.vtkdatawritesize, leave=True, unit='byte', unit_scale=True, desc='Writing snapshot file {} of {}, {}'.format(i + 1, len(G.snapshots), os.path.split(snap.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not G.progressbars)
pbar = tqdm(total=snap.vtkdatawritesize, leave=True, unit='byte', unit_scale=True, desc='Writing snapshot file {} of {}, {}'.format(i + 1, len(G.snapshots), os.path.split(snap.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not progressbars)
snap.write_vtk_imagedata(pbar, G)
pbar.close()
if G.messages: print()
if messages: print()
if G.messages:
if G.gpu is None:
print('Memory (RAM) used: ~{}'.format(human_size(p.memory_info().rss)))
else:
print('Memory (RAM) used: ~{} host + ~{} GPU'.format(human_size(p.memory_info().rss), human_size(memsolve)))
if messages:
print('Memory (RAM) used: ~{}'.format(human_size(p.memory_full_info().uss)))
if gpu:
print('GPU memory (RAM) used: ~{}'.format(human_size(memsolve)))
print('Solving time [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=tsolve)))
# If geometry information to be reused between model runs then FDTDGrid
@@ -414,12 +410,12 @@ def solve_cpu(currentmodelrun, modelend, G):
G (class): Grid class instance - holds essential parameters describing the model.
Returns:
tsolve (float): Time taken to execute solving
tsolve (float): Time taken to execute solving (seconds)
"""
tsolvestart = timer()
for iteration in tqdm(range(G.iterations), desc='Running simulation, model ' + str(currentmodelrun) + '/' + str(modelend), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not G.progressbars):
for iteration in tqdm(range(G.iterations), desc='Running simulation, model ' + str(currentmodelrun) + '/' + str(modelend), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not progressbars):
# Store field component values for every receiver and transmission line
store_outputs(iteration, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G)
@@ -429,7 +425,7 @@ def solve_cpu(currentmodelrun, modelend, G):
snap.store(G)
# Update magnetic field components
update_magnetic(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsH, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
update_magnetic(G.nx, G.ny, G.nz, hostinfo['ompthreads'], G.updatecoeffsH, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
# Update magnetic field components with the PML correction
for pml in G.pmls:
@@ -442,13 +438,13 @@ def solve_cpu(currentmodelrun, modelend, G):
# Update electric field components
# All materials are non-dispersive so do standard update
if Material.maxpoles == 0:
update_electric(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
update_electric(G.nx, G.ny, G.nz, hostinfo['ompthreads'], G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
# 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).
elif Material.maxpoles == 1:
update_electric_dispersive_1pole_A(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsE, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
update_electric_dispersive_1pole_A(G.nx, G.ny, G.nz, hostinfo['ompthreads'], G.updatecoeffsE, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
elif Material.maxpoles > 1:
update_electric_dispersive_multipole_A(G.nx, G.ny, G.nz, G.nthreads, Material.maxpoles, G.updatecoeffsE, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
update_electric_dispersive_multipole_A(G.nx, G.ny, G.nz, hostinfo['ompthreads'], Material.maxpoles, G.updatecoeffsE, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz)
# Update electric field components with the PML correction
for pml in G.pmls:
@@ -463,9 +459,9 @@ def solve_cpu(currentmodelrun, modelend, G):
# field values). Therefore it can only be completely updated after the
# electric field has been updated by the PML and source updates.
if Material.maxpoles == 1:
update_electric_dispersive_1pole_B(G.nx, G.ny, G.nz, G.nthreads, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez)
update_electric_dispersive_1pole_B(G.nx, G.ny, G.nz, hostinfo['ompthreads'], G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez)
elif Material.maxpoles > 1:
update_electric_dispersive_multipole_B(G.nx, G.ny, G.nz, G.nthreads, Material.maxpoles, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez)
update_electric_dispersive_multipole_B(G.nx, G.ny, G.nz, hostinfo['ompthreads'], Material.maxpoles, G.updatecoeffsdispersive, G.ID, G.Tx, G.Ty, G.Tz, G.Ex, G.Ey, G.Ez)
tsolve = timer() - tsolvestart
@@ -481,8 +477,8 @@ def solve_gpu(currentmodelrun, modelend, G):
G (class): Grid class instance - holds essential parameters describing the model.
Returns:
tsolve (float): Time taken to execute solving
memsolve (int): memory usage on final iteration in bytes
tsolve (float): Time taken to execute solving (seconds)
memsolve (int): memory usage on final iteration (bytes)
"""
import pycuda.driver as drv
@@ -496,7 +492,7 @@ def solve_gpu(currentmodelrun, modelend, G):
compiler_opts = None
# Create device handle and context on specifc GPU device (and make it current context)
dev = drv.Device(G.gpu.deviceID)
dev = drv.Device(gpu.deviceID)
ctx = dev.make_context()
# Electric and magnetic field updates - prepare kernels, and get kernel functions
@@ -510,8 +506,8 @@ def solve_gpu(currentmodelrun, modelend, G):
# Copy material coefficient arrays to constant memory of GPU (must be <64KB) for fields kernels
updatecoeffsE = kernels_fields.get_global('updatecoeffsE')[0]
updatecoeffsH = kernels_fields.get_global('updatecoeffsH')[0]
if G.updatecoeffsE.nbytes + G.updatecoeffsH.nbytes > G.gpu.constmem:
raise GeneralError('Too many materials in the model to fit onto constant memory of size {} on {} - {} GPU'.format(human_size(G.gpu.constmem), G.gpu.deviceID, G.gpu.name))
if G.updatecoeffsE.nbytes + G.updatecoeffsH.nbytes > gpu.constmem:
raise GeneralError('Too many materials in the model to fit onto constant memory of size {} on {} - {} GPU'.format(human_size(gpu.constmem), gpu.deviceID, gpu.name))
else:
drv.memcpy_htod(updatecoeffsE, G.updatecoeffsE)
drv.memcpy_htod(updatecoeffsH, G.updatecoeffsH)
@@ -585,7 +581,7 @@ def solve_gpu(currentmodelrun, modelend, G):
iterend = drv.Event()
iterstart.record()
for iteration in tqdm(range(G.iterations), desc='Running simulation, model ' + str(currentmodelrun) + '/' + str(modelend), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not G.progressbars):
for iteration in tqdm(range(G.iterations), desc='Running simulation, model ' + str(currentmodelrun) + '/' + str(modelend), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not progressbars):
# Get GPU memory usage on final iteration
if iteration == G.iterations - 1:
@@ -602,7 +598,7 @@ def solve_gpu(currentmodelrun, modelend, G):
# Store any snapshots
for i, snap in enumerate(G.snapshots):
if snap.time == iteration + 1:
if not G.snapsgpu2cpu:
if not snapsgpu2cpu:
store_snapshot_gpu(np.int32(i), np.int32(snap.xs),
np.int32(snap.xf), np.int32(snap.ys),
np.int32(snap.yf), np.int32(snap.zs),
@@ -698,7 +694,7 @@ def solve_gpu(currentmodelrun, modelend, G):
gpu_get_rx_array(rxs_gpu.get(), rxcoords_gpu.get(), G)
# Copy data from any snapshots back to correct snapshot objects
if G.snapshots and not G.snapsgpu2cpu:
if G.snapshots and not snapsgpu2cpu:
for i, snap in enumerate(G.snapshots):
gpu_get_snapshot_array(snapEx_gpu.get(), snapEy_gpu.get(), snapEz_gpu.get(),
snapHx_gpu.get(), snapHy_gpu.get(), snapHz_gpu.get(), i, snap)

490
gprMax/optimisation_taguchi.py
查看文件

@@ -1,490 +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 <http://www.gnu.org/licenses/>.
from collections import OrderedDict
import datetime
from importlib import import_module
import os
import pickle
import sys
from time import perf_counter
from colorama import init, Fore, Style
init()
import numpy as np
from gprMax.constants import floattype
from gprMax.exceptions import CmdInputError
from gprMax.gprMax import run_std_sim
from gprMax.gprMax import run_mpi_sim
from gprMax.utilities import get_terminal_width
from gprMax.utilities import open_path_file
def run_opt_sim(args, inputfile, usernamespace):
"""Run a simulation using Taguchi's optmisation process.
Args:
args (dict): Namespace with command line arguments
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.
"""
tsimstart = perf_counter()
if args.n > 1:
raise CmdInputError('When a Taguchi optimisation is being carried out the number of model runs argument is not required')
inputfileparts = os.path.splitext(inputfile.name)
# Default maximum number of iterations of optimisation to perform (used
# if the stopping criterion is not achieved)
maxiterations = 20
# Process Taguchi code blocks in the input file; pass in ordered
# dictionary to hold parameters to optimise
tmp = usernamespace.copy()
tmp.update({'optparams': OrderedDict()})
taguchinamespace = taguchi_code_blocks(inputfile, tmp)
# Extract dictionaries and variables containing initialisation parameters
optparams = taguchinamespace['optparams']
fitness = taguchinamespace['fitness']
if 'maxiterations' in taguchinamespace:
maxiterations = taguchinamespace['maxiterations']
# Store initial parameter ranges
optparamsinit = list(optparams.items())
# Dictionary to hold history of optmised values of parameters
optparamshist = OrderedDict((key, list()) for key in optparams)
# Import specified fitness function
fitness_metric = getattr(import_module('user_libs.optimisation_taguchi.fitness_functions'), fitness['name'])
# Select OA
OA, N, cols, k, s, t = construct_OA(optparams)
taguchistr = '\n--- Taguchi optimisation'
print('{} {}\n'.format(taguchistr, '-' * (get_terminal_width() - 1 - len(taguchistr))))
print('Orthogonal array: {:g} experiments per iteration, {:g} parameters ({:g} will be used), {:g} levels, and strength {:g}'.format(N, cols, k, s, t))
tmp = [(k, v) for k, v in optparams.items()]
print('Parameters to optimise with ranges: {}'.format(str(tmp).strip('[]')))
print('Output name(s) from model: {}'.format(fitness['args']['outputs']))
print('Fitness function "{}" with stopping criterion {:g}'.format(fitness['name'], fitness['stop']))
print('Maximum iterations: {:g}'.format(maxiterations))
# Initialise arrays and lists to store parameters required throughout optimisation
# Lower, central, and upper values for each parameter
levels = np.zeros((s, k), dtype=floattype)
# Optimal lower, central, or upper value for each parameter
levelsopt = np.zeros(k, dtype=np.uint8)
# Difference used to set values for levels
levelsdiff = np.zeros(k, dtype=floattype)
# History of fitness values from each confirmation experiment
fitnessvalueshist = []
iteration = 0
while iteration < maxiterations:
# Reset number of model runs to number of experiments
args.n = N
usernamespace['number_model_runs'] = N
# Fitness values for each experiment
fitnessvalues = []
# Set parameter ranges and define experiments
optparams, levels, levelsdiff = calculate_ranges_experiments(optparams, optparamsinit, levels, levelsopt, levelsdiff, OA, N, k, s, iteration)
# Run model for each experiment
# Mixed mode MPI with OpenMP or CUDA - MPI task farm for models with
# each model parallelised with OpenMP (CPU) or CUDA (GPU)
if args.mpi:
run_mpi_sim(args, inputfile, usernamespace, optparams)
# Standard behaviour - models run serially with each model parallelised
# with OpenMP (CPU) or CUDA (GPU)
else:
run_std_sim(args, inputfile, usernamespace, optparams)
# Calculate fitness value for each experiment
for experiment in range(1, N + 1):
outputfile = inputfileparts[0] + str(experiment) + '.out'
fitnessvalues.append(fitness_metric(outputfile, fitness['args']))
os.remove(outputfile)
taguchistr = '\n--- Taguchi optimisation, iteration {}: {} initial experiments with fitness values {}.'.format(iteration + 1, N, fitnessvalues)
print('{} {}\n'.format(taguchistr, '-' * (get_terminal_width() - 1 - len(taguchistr))))
# Calculate optimal levels from fitness values by building a response
# table; update dictionary of parameters with optimal values
optparams, levelsopt = calculate_optimal_levels(optparams, levels, levelsopt, fitnessvalues, OA, N, k)
# Update dictionary with history of parameters with optimal values
for key, value in optparams.items():
optparamshist[key].append(value[0])
# Run a confirmation experiment with optimal values
args.n = 1
usernamespace['number_model_runs'] = 1
# Mixed mode MPI with OpenMP or CUDA - MPI task farm for models with
# each model parallelised with OpenMP (CPU) or CUDA (GPU)
if args.mpi:
run_mpi_sim(args, inputfile, usernamespace, optparams)
# Standard behaviour - models run serially with each model parallelised
# with OpenMP (CPU) or CUDA (GPU)
else:
run_std_sim(args, inputfile, usernamespace, optparams)
# Calculate fitness value for confirmation experiment
outputfile = inputfileparts[0] + '.out'
fitnessvalueshist.append(fitness_metric(outputfile, fitness['args']))
# Rename confirmation experiment output file so that it is retained for each iteraction
os.rename(outputfile, os.path.splitext(outputfile)[0] + '_final' + str(iteration + 1) + '.out')
taguchistr = '\n--- Taguchi optimisation, iteration {} completed. History of optimal parameter values {} and of fitness values {}'.format(iteration + 1, dict(optparamshist), fitnessvalueshist)
print('{} {}\n'.format(taguchistr, '-' * (get_terminal_width() - 1 - len(taguchistr))))
iteration += 1
# Stop optimisation if stopping criterion has been reached
if fitnessvalueshist[iteration - 1] > fitness['stop']:
taguchistr = '\n--- Taguchi optimisation stopped as fitness criteria reached: {:g} > {:g}'.format(fitnessvalueshist[iteration - 1], fitness['stop'])
print('{} {}\n'.format(taguchistr, '-' * (get_terminal_width() - 1 - len(taguchistr))))
break
# Stop optimisation if successive fitness values are within a percentage threshold
fitnessvaluesthres = 0.1
if iteration > 2:
fitnessvaluesclose = (np.abs(fitnessvalueshist[iteration - 2] - fitnessvalueshist[iteration - 1]) / fitnessvalueshist[iteration - 1]) * 100
if fitnessvaluesclose < fitnessvaluesthres:
taguchistr = '\n--- Taguchi optimisation stopped as successive fitness values within {}%'.format(fitnessvaluesthres)
print('{} {}\n'.format(taguchistr, '-' * (get_terminal_width() - 1 - len(taguchistr))))
break
tsimend = perf_counter()
# Save optimisation parameters history and fitness values history to file
opthistfile = inputfileparts[0] + '_hist.pickle'
with open(opthistfile, 'wb') as f:
pickle.dump(optparamshist, f)
pickle.dump(fitnessvalueshist, f)
pickle.dump(optparamsinit, f)
taguchistr = '\n=== Taguchi optimisation completed in [HH:MM:SS]: {} after {} iteration(s)'.format(datetime.timedelta(seconds=int(tsimend - tsimstart)), iteration)
print('{} {}\n'.format(taguchistr, '=' * (get_terminal_width() - 1 - len(taguchistr))))
print('History of optimal parameter values {} and of fitness values {}\n'.format(dict(optparamshist), fitnessvalueshist))
def taguchi_code_blocks(inputfile, taguchinamespace):
"""
Looks for and processes a Taguchi code block (containing Python code) in
the input file. It will ignore any lines that are comments, i.e. begin
with a double hash (##), and any blank lines.
Args:
inputfile (object): File object for the input file.
taguchinamespace (dict): Namespace that can be accessed by user a
Taguchi code block in input file.
Returns:
processedlines (list): Input commands after Python processing.
"""
# Strip out any newline characters and comments that must begin with double hashes
inputlines = [line.rstrip() for line in inputfile if(not line.startswith('##') and line.rstrip('\n'))]
# Rewind input file in preparation for passing to standard command reading function
inputfile.seek(0)
# Store length of dict
taglength = len(taguchinamespace)
x = 0
while(x < len(inputlines)):
if(inputlines[x].startswith('#taguchi:')):
# String to hold Python code to be executed
taguchicode = ''
x += 1
while not inputlines[x].startswith('#end_taguchi:'):
# Add all code in current code block to string
taguchicode += inputlines[x] + '\n'
x += 1
if x == len(inputlines):
raise CmdInputError('Cannot find the end of the Taguchi code block, i.e. missing #end_taguchi: command.')
# Compile code for faster execution
taguchicompiledcode = compile(taguchicode, '<string>', 'exec')
# Execute code block & make available only usernamespace
exec(taguchicompiledcode, taguchinamespace)
x += 1
# Check if any Taguchi code blocks were found
if len(taguchinamespace) == taglength:
raise CmdInputError('No #taguchi and #end_taguchi code blocks found.')
return taguchinamespace
def construct_OA(optparams):
"""
Load an orthogonal array (OA) from a numpy file. Configure and
return OA and properties of OA.
Args:
optparams (dict): Dictionary containing name of parameters to
optimise and their initial ranges
Returns:
OA (array): Orthogonal array
N (int): Number of experiments in OA
cols (int): Number of columns in OA
k (int): Number of columns in OA cut down to number of parameters to optimise
s (int): Number of levels in OA
t (int): Strength of OA
"""
oadirectory = os.path.join(os.path.dirname(os.path.abspath(__file__)), os.pardir, 'user_libs', 'optimisation_taguchi')
oadirectory = os.path.abspath(oadirectory)
# Properties of the orthogonal array (OA)
# Strength
t = 2
# Number of levels
s = 3
# Number of parameters to optimise
k = len(optparams)
# Load the appropriate OA
if k <= 4:
OA = np.load(os.path.join(oadirectory, 'OA_9_4_3_2.npy'))
# Number of experiments
N = OA.shape[0]
# Number of columns of OA before cut down
cols = OA.shape[1]
# Cut down OA columns to number of parameters to optimise
OA = OA[:, 0:k]
elif k <= 7:
OA = np.load(os.path.join(oadirectory, 'OA_18_7_3_2.npy'))
# Number of experiments
N = OA.shape[0]
# Number of columns of OA before cut down
cols = OA.shape[1]
# Cut down OA columns to number of parameters to optimise
OA = OA[:, 0:k]
else:
# THIS CASE NEEDS FURTHER TESTING
print(Fore.RED + 'WARNING: Optimising more than 7 parameters is currently an experimental feature!' + Style.RESET_ALL)
p = int(np.ceil(np.log(k * (s - 1) + 1) / np.log(s)))
# Number of experiments
N = s**p
# Number of columns
cols = int((N - 1) / (s - 1))
# Algorithm to construct OA from:
# http://ieeexplore.ieee.org/xpl/articleDetails.jsp?reload=true&arnumber=6812898
OA = np.zeros((N + 1, cols + 1), dtype=np.int8)
# Construct basic columns
for ii in range(1, p + 1):
col = int((s**(ii - 1) - 1) / (s - 1) + 1)
for row in range(1, N + 1):
OA[row, col] = np.mod(np.floor((row - 1) / (s**(p - ii))), s)
# Construct non-basic columns
for ii in range(2, p + 1):
col = int((s**(ii - 1) - 1) / (s - 1) + 1)
for jj in range(1, col):
for kk in range(1, s):
OA[:, col + (jj - 1) * (s - 1) + kk] = np.mod(OA[:, jj] * kk + OA[:, col], s)
# First row and first columns are unneccessary, only there to
# match algorithm, and cut down columns to number of parameters to optimise
OA = OA[1:, 1:k + 1]
return OA, N, cols, k, s, t
def calculate_ranges_experiments(optparams, optparamsinit, levels, levelsopt, levelsdiff, OA, N, k, s, i):
"""Calculate values for parameters to optimise for a set of experiments.
Args:
optparams (dict): Ordered dictionary containing name of parameters to optimise and their values
optparamsinit (list): Initial ranges for parameters to optimise
levels (array): Lower, central, and upper values for each parameter
levelsopt (array): Optimal level for each parameter from previous iteration
levelsdiff (array): Difference used to set values in levels array
OA (array): Orthogonal array
N (int): Number of experiments in OA
k (int): Number of parameters to optimise in OA
s (int): Number of levels in OA
i (int): Iteration number
Returns:
optparams (dict): Ordered dictionary containing name of parameters to optimise and their values
levels (array): Lower, central, and upper values for each parameter
levelsdiff (array): Difference used to set values in levels array
"""
# Gaussian reduction function used for calculating levels
T = 18 # Usually values between 15 - 20
RR = np.exp(-(i / T)**2)
# Calculate levels for each parameter
for p in range(k):
# Set central level for first iteration to midpoint of initial range and don't use RR
if i == 0:
levels[1, p] = ((optparamsinit[p][1][1] - optparamsinit[p][1][0]) / 2) + optparamsinit[p][1][0]
levelsdiff[p] = (optparamsinit[p][1][1] - optparamsinit[p][1][0]) / (s + 1)
# Set central level to optimum from previous iteration
else:
levels[1, p] = levels[levelsopt[p], p]
levelsdiff[p] = RR * levelsdiff[p]
# Set levels if below initial range
if levels[1, p] - levelsdiff[p] < optparamsinit[p][1][0]:
levels[0, p] = optparamsinit[p][1][0]
levels[1, p] = optparamsinit[p][1][0] + levelsdiff[p]
levels[2, p] = optparamsinit[p][1][0] + 2 * levelsdiff[p]
# Set levels if above initial range
elif levels[1, p] + levelsdiff[p] > optparamsinit[p][1][1]:
levels[0, p] = optparamsinit[p][1][1] - 2 * levelsdiff[p]
levels[1, p] = optparamsinit[p][1][1] - levelsdiff[p]
levels[2, p] = optparamsinit[p][1][1]
# Set levels normally
else:
levels[0, p] = levels[1, p] - levelsdiff[p]
levels[2, p] = levels[1, p] + levelsdiff[p]
# Update dictionary of parameters to optimise with lists of new values; clear dictionary first
optparams = OrderedDict((key, list()) for key in optparams)
p = 0
for key, value in optparams.items():
for exp in range(N):
if OA[exp, p] == 0:
optparams[key].append(levels[0, p])
elif OA[exp, p] == 1:
optparams[key].append(levels[1, p])
elif OA[exp, p] == 2:
optparams[key].append(levels[2, p])
p += 1
return optparams, levels, levelsdiff
def calculate_optimal_levels(optparams, levels, levelsopt, fitnessvalues, OA, N, k):
"""Calculate optimal levels from results of fitness metric by building a response table.
Args:
optparams (dict): Ordered dictionary containing name of parameters to optimise and their values
levels (array): Lower, central, and upper values for each parameter
levelsopt (array): Optimal level for each parameter from previous iteration
fitnessvalues (list): Values from results of fitness metric
OA (array): Orthogonal array
N (int): Number of experiments in OA
k (int): Number of parameters to optimise in OA
Returns:
optparams (dict): Ordered dictionary containing name of parameters to optimise and their values
levelsopt (array): Optimal level for each parameter from previous iteration
"""
# Build a table of responses based on the results of the fitness metric
for p in range(k):
responses = np.zeros(3, dtype=floattype)
cnt1 = 0
cnt2 = 0
cnt3 = 0
for exp in range(N):
if OA[exp, p] == 0:
responses[0] += fitnessvalues[exp]
cnt1 += 1
elif OA[exp, p] == 1:
responses[1] += fitnessvalues[exp]
cnt2 += 1
elif OA[exp, p] == 2:
responses[2] += fitnessvalues[exp]
cnt3 += 1
responses[0] /= cnt1
responses[1] /= cnt2
responses[2] /= cnt3
# Calculate optimal level from table of responses
optlevel = np.where(responses == np.amax(responses))[0]
# If 2 experiments produce the same fitness value pick first level
# (this shouldn't happen if the fitness function is designed correctly)
if len(optlevel) > 1:
optlevel = optlevel[0]
levelsopt[p] = optlevel
# Update dictionary of parameters to optimise with lists of new values; clear dictionary first
optparams = OrderedDict((key, list()) for key in optparams)
p = 0
for key, value in optparams.items():
optparams[key].append(levels[levelsopt[p], p])
p += 1
return optparams, levelsopt
def plot_optimisation_history(fitnessvalueshist, optparamshist, optparamsinit):
"""Plot the history of fitness values and each optimised parameter values for the optimisation.
Args:
fitnessvalueshist (list): History of fitness values
optparamshist (dict): Name of parameters to optimise and history of their values
"""
import matplotlib.pyplot as plt
# Plot history of fitness values
fig, ax = plt.subplots(subplot_kw=dict(xlabel='Iterations', ylabel='Fitness value'), num='History of fitness values', figsize=(20, 10), facecolor='w', edgecolor='w')
iterations = np.arange(1, len(fitnessvalueshist) + 1)
ax.plot(iterations, fitnessvalueshist, 'r', marker='.', ms=15, lw=1)
ax.set_xlim(1, len(fitnessvalueshist))
ax.grid()
# Plot history of optimisation parameters
p = 0
for key, value in optparamshist.items():
fig, ax = plt.subplots(subplot_kw=dict(xlabel='Iterations', ylabel='Parameter value'), num='History of ' + key + ' parameter', figsize=(20, 10), facecolor='w', edgecolor='w')
ax.plot(iterations, optparamshist[key], 'r', marker='.', ms=15, lw=1)
ax.set_xlim(1, len(fitnessvalueshist))
ax.set_ylim(optparamsinit[p][1][0], optparamsinit[p][1][1])
ax.grid()
p += 1
plt.show()

11
gprMax/pml.py
查看文件

@@ -21,9 +21,10 @@ from importlib import import_module
import numpy as np
from tqdm import tqdm
from gprMax.constants import e0
from gprMax.constants import z0
from gprMax.constants import floattype
from gprMax.config import e0
from gprMax.config import z0
from gprMax.config import floattype
from gprMax.config import hostinfo
from gprMax.exceptions import GeneralError
@@ -281,7 +282,7 @@ class PML(object):
pmlmodule = 'gprMax.pml_updates.pml_updates_electric_' + G.pmlformulation + '_ext'
func = getattr(import_module(pmlmodule), 'order' + str(len(self.CFS)) + '_' + self.direction)
func(self.xs, self.xf, self.ys, self.yf, self.zs, self.zf, G.nthreads, G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, self.EPhi1, self.EPhi2, self.ERA, self.ERB, self.ERE, self.ERF, self.d)
func(self.xs, self.xf, self.ys, self.yf, self.zs, self.zf, hostinfo['ompthreads'], G.updatecoeffsE, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, self.EPhi1, self.EPhi2, self.ERA, self.ERB, self.ERE, self.ERF, self.d)
def update_magnetic(self, G):
"""This functions updates magnetic field components with the PML correction.
@@ -292,7 +293,7 @@ class PML(object):
pmlmodule = 'gprMax.pml_updates.pml_updates_magnetic_' + G.pmlformulation + '_ext'
func = getattr(import_module(pmlmodule), 'order' + str(len(self.CFS)) + '_' + self.direction)
func(self.xs, self.xf, self.ys, self.yf, self.zs, self.zf, G.nthreads, G.updatecoeffsH, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, self.HPhi1, self.HPhi2, self.HRA, self.HRB, self.HRE, self.HRF, self.d)
func(self.xs, self.xf, self.ys, self.yf, self.zs, self.zf, hostinfo['ompthreads'], G.updatecoeffsH, G.ID, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, self.HPhi1, self.HPhi2, self.HRA, self.HRB, self.HRE, self.HRF, self.d)
def gpu_set_blocks_per_grid(self, G):
"""Set the blocks per grid size used for updating the PML field arrays on a GPU.

2
gprMax/pml_updates/pml_updates_electric_HORIPML_ext.pyx
查看文件

@@ -20,7 +20,7 @@ import numpy as np
cimport numpy as np
from cython.parallel import prange
from gprMax.constants cimport floattype_t
from gprMax.config cimport floattype_t
cpdef void order1_xminus(

2
gprMax/pml_updates/pml_updates_electric_MRIPML_ext.pyx
查看文件

@@ -20,7 +20,7 @@ import numpy as np
cimport numpy as np
from cython.parallel import prange
from gprMax.constants cimport floattype_t
from gprMax.config cimport floattype_t
cpdef void order1_xminus(

2
gprMax/pml_updates/pml_updates_magnetic_HORIPML_ext.pyx
查看文件

@@ -20,7 +20,7 @@ import numpy as np
cimport numpy as np
from cython.parallel import prange
from gprMax.constants cimport floattype_t
from gprMax.config cimport floattype_t
cpdef void order1_xminus(

2
gprMax/pml_updates/pml_updates_magnetic_MRIPML_ext.pyx
查看文件

@@ -20,7 +20,7 @@ import numpy as np
cimport numpy as np
from cython.parallel import prange
from gprMax.constants cimport floattype_t
from gprMax.config cimport floattype_t
cpdef void order1_xminus(

2
gprMax/receivers.py
查看文件

@@ -20,7 +20,7 @@ from collections import OrderedDict
import numpy as np
from gprMax.constants import floattype
from gprMax.config import floattype
class Rx(object):

5
gprMax/snapshots.py
查看文件

@@ -22,7 +22,8 @@ from struct import pack
import numpy as np
from gprMax.constants import floattype
from gprMax.config import floattype
from gprMax.config import snapsgpu2cpu
from gprMax.snapshots_ext import calculate_snapshot_fields
from gprMax.utilities import round_value
@@ -192,7 +193,7 @@ def gpu_initialise_snapshot_array(G):
Snapshot.bpg = (int(np.ceil(((Snapshot.nx_max) * (Snapshot.ny_max) * (Snapshot.nz_max)) / Snapshot.tpb[0])), 1, 1)
# 4D arrays to store snapshots on GPU, e.g. snapEx(time, x, y, z)
numsnaps = 1 if G.snapsgpu2cpu else len(G.snapshots)
numsnaps = 1 if snapsgpu2cpu else len(G.snapshots)
snapEx = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=floattype)
snapEy = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=floattype)
snapEz = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max), dtype=floattype)

2
gprMax/snapshots_ext.pyx
查看文件

@@ -18,7 +18,7 @@
cimport numpy as np
from gprMax.constants cimport floattype_t
from gprMax.config cimport floattype_t
cpdef void calculate_snapshot_fields(

4
gprMax/sources.py
查看文件

@@ -20,8 +20,8 @@ from copy import deepcopy
import numpy as np
from gprMax.constants import c
from gprMax.constants import floattype
from gprMax.config import c
from gprMax.config import floattype
from gprMax.grid import Ix
from gprMax.grid import Iy
from gprMax.grid import Iz

8
gprMax/utilities.py
查看文件

@@ -35,10 +35,7 @@ init()
import numpy as np
from time import process_time
from gprMax.constants import complextype
from gprMax.constants import floattype
from gprMax.exceptions import GeneralError
from gprMax.materials import Material
def get_terminal_width():
@@ -335,6 +332,7 @@ def get_host_info():
if not hostinfo['physicalcores']:
hostinfo['physicalcores'] = hostinfo['logicalcores']
hostinfo['ram'] = psutil.virtual_memory().total
hostinfo['ompthreads'] = 1
return hostinfo
@@ -353,6 +351,10 @@ class GPU(object):
self.pcibusID = None
self.constmem = None
self.totalmem = None
# Threads per block for main field updates
self.tpb = (256, 1, 1)
# Blocks per grid for main field updates (set in grid.py)
self.bpg = None
def get_gpu_info(self, drv):
"""Set information about GPU.

3
tests/analytical_solutions.py
查看文件

@@ -1,6 +1,7 @@
import numpy as np
from gprMax.constants import c, e0
from gprMax.config import c
from gprMax.config import e0
from gprMax.waveforms import Waveform

3
user_libs/antenna_patterns/initial_save.py
查看文件

@@ -13,7 +13,8 @@ import h5py
import numpy as np
import matplotlib.pyplot as plt
from gprMax.constants import c, z0
from gprMax.config import c
from gprMax.config import z0
# Parse command line arguments

5
user_libs/antenna_patterns/plot_fields.py
查看文件

@@ -12,7 +12,8 @@ import sys
import numpy as np
import matplotlib.pyplot as plt
from gprMax.constants import c, z0
from gprMax.config import c
from gprMax.config import z0
# Parse command line arguments
@@ -89,7 +90,7 @@ for patt in range(0, len(radii)):
power = 10 * np.log10(pattplot)
# Replace any NaNs or Infs from zero division
power[np.invert(np.isfinite(power))] = 0
ax.plot(theta, power, label='{:.2f}m'.format(radii[patt]), marker='.', ms=6, lw=1.5)
# Add Hertzian dipole plot

二进制
user_libs/optimisation_taguchi/OA_18_7_3_2.npy
查看文件

二进制文件未显示。

二进制
user_libs/optimisation_taguchi/OA_9_4_3_2.npy
查看文件

二进制文件未显示。

92
user_libs/optimisation_taguchi/antenna_bowtie_opt.in
查看文件

@@ -1,92 +0,0 @@
#taguchi:
optparams['resinner'] = [0.1, 1000]
optparams['resmiddle'] = [0.1, 1000]
optparams['resouter'] = [0.1, 1000]
fitness = {'name': 'min_max_value', 'stop': 10, 'args': {'type': 'absmax', 'outputs': 'Ex60mm'}}
#end_taguchi:
#python:
import numpy as np
from gprMax.input_cmd_funcs import *
title = 'antenna_bowtie_opt'
print('#title: {}'.format(title))
domain = domain(0.180, 0.120, 0.160)
dxdydz = dx_dy_dz(0.001, 0.001, 0.001)
timewindow = time_window(5e-9)
fr4_dims = (0.120, 0.060, 0.002)
bowtie_dims = (0.050, 0.040) # Length, height
flare_angle = np.arctan((bowtie_dims[1]/2) / bowtie_dims[0])
tx_pos = (domain[0]/2, domain[1]/2, 0.050)
# Vertical slot positions, relative to feed position, i.e. txpos[0]
vcut_pos = (0.014, 0.027, 0.038)
# Loading resistor values
res = np.array([optparams['resinner'], optparams['resmiddle'], optparams['resouter']])
rescond = ((1 / res) * (dxdydz[1] / (dxdydz[0] * dxdydz[2]))) / 2 # Divide by number of parallel edges per resistor
# Materials
material(4.8, 0, 1, 0, 'fr4')
for i in range(len(res)):
material(1, rescond[i], 1, 0, 'res' + str(i + 1))
# Source excitation and type
print('#waveform: gaussian 1 2e9 mypulse')
print('#transmission_line: x {:g} {:g} {:g} 50 mypulse'.format(tx_pos[0], tx_pos[1], tx_pos[2]))
# Output point - distance from tx_pos in z direction
print('#rx: {:g} {:g} {:g} Ex60mm Ex'.format(tx_pos[0], tx_pos[1], tx_pos[2] + 0.060))
# Bowtie - upper x half
triangle(tx_pos[0], tx_pos[1], tx_pos[2], tx_pos[0] + bowtie_dims[0], tx_pos[1] - bowtie_dims[1]/2, tx_pos[2], tx_pos[0] + bowtie_dims[0], tx_pos[1] + bowtie_dims[1]/2, tx_pos[2], 0, 'pec')
# Bowtie - upper x half - vertical cuts
for i in range(len(vcut_pos)):
for j in range(int(bowtie_dims[1] / dxdydz[2])):
edge(tx_pos[0] + vcut_pos[i], tx_pos[1] - bowtie_dims[1]/2 + j * dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] - bowtie_dims[1]/2 + j * dxdydz[1], tx_pos[2], 'free_space')
# Bowtie - upper x half - vertical cuts - loading
for i in range(len(vcut_pos)):
gap = ((vcut_pos[i] * np.tan(flare_angle) * 2) - 4*dxdydz[1]) / 5
edge(tx_pos[0] + vcut_pos[i], tx_pos[1] - (1.5 * gap) - dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] - (1.5 * gap) - dxdydz[1], tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] + vcut_pos[i], tx_pos[1] - (1.5 * gap) - 2*dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] - (1.5 * gap) - 2*dxdydz[1], tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] + vcut_pos[i], tx_pos[1] - (0.5 * gap), tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] - (0.5 * gap), tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] + vcut_pos[i], tx_pos[1] - (0.5 * gap) - dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] - (0.5 * gap) - dxdydz[1], tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] + vcut_pos[i], tx_pos[1] + (0.5 * gap), tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] + (0.5 * gap), tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] + vcut_pos[i], tx_pos[1] + (0.5 * gap) + dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] + (0.5 * gap) + dxdydz[1], tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] + vcut_pos[i], tx_pos[1] + (1.5 * gap) + dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] + (1.5 * gap) + dxdydz[1], tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] + vcut_pos[i], tx_pos[1] + (1.5 * gap) + 2*dxdydz[1], tx_pos[2], tx_pos[0] + vcut_pos[i] + dxdydz[0], tx_pos[1] + (1.5 * gap) + 2*dxdydz[1], tx_pos[2], 'res' + str(i + 1))
# Bowtie - lower x half
triangle(tx_pos[0] + dxdydz[0], tx_pos[1], tx_pos[2], tx_pos[0] - bowtie_dims[0], tx_pos[1] - bowtie_dims[1]/2, tx_pos[2], tx_pos[0] - bowtie_dims[0], tx_pos[1] + bowtie_dims[1]/2, tx_pos[2], 0, 'pec')
# Bowtie - lower x half - cuts for loading
for i in range(len(vcut_pos)):
for j in range(int(bowtie_dims[1] / dxdydz[2])):
edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] - bowtie_dims[1]/2 + j * dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] - bowtie_dims[1]/2 + j * dxdydz[1], tx_pos[2], 'free_space')
# Bowtie - lower x half - vertical cuts - loading
for i in range(len(vcut_pos)):
gap = ((vcut_pos[i] * np.tan(flare_angle) * 2) - 4*dxdydz[1]) / 5
edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] - (1.5 * gap) - dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] - (1.5 * gap) - dxdydz[1], tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] - (1.5 * gap) - 2*dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] - (1.5 * gap) - 2*dxdydz[1], tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] - (0.5 * gap), tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] - (0.5 * gap), tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] - (0.5 * gap) - dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] - (0.5 * gap) - dxdydz[1], tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] + (0.5 * gap), tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] + (0.5 * gap), tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] + (0.5 * gap) + dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] + (0.5 * gap) + dxdydz[1], tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] + (1.5 * gap) + dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] + (1.5 * gap) + dxdydz[1], tx_pos[2], 'res' + str(i + 1))
edge(tx_pos[0] - vcut_pos[i] - dxdydz[0], tx_pos[1] + (1.5 * gap) + 2*dxdydz[1], tx_pos[2], tx_pos[0] - vcut_pos[i], tx_pos[1] + (1.5 * gap) + 2*dxdydz[1], tx_pos[2], 'res' + str(i + 1))
# PCB
box(tx_pos[0] - fr4_dims[0]/2, tx_pos[1] - fr4_dims[1]/2, tx_pos[2] - fr4_dims[2], tx_pos[0] + fr4_dims[0]/2, tx_pos[1] + fr4_dims[1]/2, tx_pos[2], 'fr4')
# Detailed geometry view of PCB and bowtie
#geometry_view(tx_pos[0] - fr4_dims[0]/2, tx_pos[1] - fr4_dims[1]/2, tx_pos[2] - fr4_dims[2], tx_pos[0] + fr4_dims[0]/2, tx_pos[1] + fr4_dims[1]/2, tx_pos[2] + dxdydz[2], dxdydz[0], dxdydz[1], dxdydz[2], title + '_tx', type='f')
# Geometry view of entire domain
#geometry_view(0, 0, 0, domain[0], domain[1], domain[2], dxdydz[0], dxdydz[1], dxdydz[2], title)
#end_python:

362
user_libs/optimisation_taguchi/fitness_functions.py
查看文件

@@ -1,362 +0,0 @@
# Copyright (C) 2015-2019, Craig Warren
#
# This module is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License.
# To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/.
#
# Please use the attribution at http://dx.doi.org/10.1190/1.3548506
import os
import sys
import h5py
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from gprMax.exceptions import GeneralError
"""This module contains fitness metric functions that can be used with the Taguchi optimisation method.
All fitness functions must take two arguments and return a single fitness value.
The first argument should be the name of the output file
The second argument is a dictionary which can contain any number of additional arguments, e.g. names (IDs) of outputs (rxs) from input file
"""
def min_max_value(filename, args):
"""Minimum value from a response.
Args:
filename (str): Name of output file
args (dict): 'type' key with string 'min', 'max' or 'absmax'; 'outputs' key with a list of names (IDs) of outputs (rxs) from input file
Returns:
value (float): Minimum, maximum, or absolute maximum value from specific outputs
"""
f = h5py.File(filename, 'r')
nrx = f.attrs['nrx']
value = 0
outputsused = False
for rx in range(1, nrx + 1):
output = f['/rxs/rx' + str(rx) + '/']
if output.attrs['Name'] in args['outputs']:
outputname = list(output.keys())[0]
if args['type'] == 'min':
value += np.abs(np.amin(output[outputname]))
elif args['type'] == 'max':
value += np.amax(output[outputname])
elif args['type'] == 'absmax':
value += np.amax(np.abs(output[outputname]))
else:
raise GeneralError('type must be either min, max, or absmax')
outputsused = True
# Check in case no outputs where found
if not outputsused:
raise GeneralError('No outputs matching {} were found'.format(args['outputs']))
return value
def xcorr(filename, args):
"""Maximum value of a cross-correlation between a response and a reference response.
Args:
filename (str): Name of output file
args (dict): 'refresp' key with path & filename of reference response (time, amp) stored in a text file; 'outputs' key with a list of names (IDs) of outputs (rxs) from input file
Returns:
xcorrmax (float): Maximum value from specific outputs
"""
# Load (from text file) the reference response. See if file exists at specified path and if not try input file directory
refrespfile = os.path.abspath(args['refresp'])
if not os.path.isfile(refrespfile):
raise GeneralError('Cannot load reference response from {}'.format(refrespfile))
with open(refrespfile, 'r') as f:
refdata = np.loadtxt(f)
reftime = refdata[:, 0] * 1e-9
refresp = refdata[:, 1]
# Load response from output file
f = h5py.File(filename, 'r')
nrx = f.attrs['nrx']
modeltime = np.arange(0, f.attrs['dt'] * f.attrs['Iterations'], f.attrs['dt'])
outputsused = False
for rx in range(1, nrx + 1):
output = f['/rxs/rx' + str(rx) + '/']
if output.attrs['Name'] in args['outputs']:
outputname = list(output.keys())[0]
modelresp = output[outputname]
# Convert electric field value (V/m) to voltage (V)
if outputname == 'Ex':
modelresp *= -f.attrs['dx, dy, dz'][0]
elif outputname == 'Ey':
modelresp *= -f.attrs['dx, dy, dz'][1]
elif outputname == 'Ez':
modelresp *= -f.attrs['dx, dy, dz'][2]
outputsused = True
# Check in case no outputs where found
if not outputsused:
raise GeneralError('No outputs matching {} were found'.format(args['outputs']))
# Normalise reference respose and response from output file
# refresp /= np.amax(np.abs(refresp))
# modelresp /= np.amax(np.abs(modelresp))
# Make both responses the same length in time
# if reftime[-1] > modeltime[-1]:
# reftime = np.arange(0, f.attrs['dt'] * f.attrs['Iterations'], reftime[-1] / len(reftime))
# refresp = refresp[0:len(reftime)]
# elif modeltime[-1] > reftime[-1]:
# modeltime = np.arange(0, reftime[-1], f.attrs['dt'])
# modelresp = modelresp[0:len(modeltime)]
# Downsample the response with the higher sampling rate
# if len(modeltime) < len(reftime):
# refresp = signal.resample(refresp, len(modelresp))
# elif len(reftime) < len(modeltime):
# modelresp = signal.resample(modelresp, len(refresp))
# Prepare data for normalized cross-correlation
refresp = (refresp - np.mean(refresp)) / (np.std(refresp) * len(refresp))
modelresp = (modelresp - np.mean(modelresp)) / np.std(modelresp)
# Plots responses for checking
# fig, ax = plt.subplots(subplot_kw=dict(xlabel='Iterations', ylabel='Voltage [V]'), figsize=(20, 10), facecolor='w', edgecolor='w')
# ax.plot(refresp,'r', lw=2, label='refresp')
# ax.plot(modelresp,'b', lw=2, label='modelresp')
# ax.grid()
# plt.show()
# Calculate cross-correlation
xcorr = signal.correlate(refresp, modelresp)
# Set any NaNs to zero
xcorr = np.nan_to_num(xcorr)
# Plot cross-correlation for checking
# fig, ax = plt.subplots(subplot_kw=dict(xlabel='Iterations', ylabel='Voltage [V]'), figsize=(20, 10), facecolor='w', edgecolor='w')
# ax.plot(xcorr,'r', lw=2, label='xcorr')
# ax.grid()
# plt.show()
xcorrmax = np.amax(xcorr)
return xcorrmax
def min_sum_diffs(filename, args):
"""Sum of the differences (in dB) between responses and a reference response.
Args:
filename (str): Name of output file
args (dict): 'refresp' key with path & filename of reference response; 'outputs' key with a list of names (IDs) of outputs (rxs) from input file
Returns:
diffdB (float): Sum of the differences (in dB) between responses and a reference response
"""
# Load (from gprMax output file) the reference response. See if file exists at specified path and if not try input file directory
refrespfile = os.path.abspath(args['refresp'])
if not os.path.isfile(refrespfile):
raise GeneralError('Cannot load reference response at {}'.format(refrespfile))
f = h5py.File(refrespfile, 'r')
tmp = f['/rxs/rx1/']
fieldname = list(tmp.keys())[0]
refresp = np.array(tmp[fieldname])
# Load (from gprMax output file) the response
f = h5py.File(filename, 'r')
nrx = f.attrs['nrx']
diffdB = 0
outputs = 0
for rx in range(1, nrx + 1):
output = f['/rxs/rx' + str(rx) + '/']
if output.attrs['Name'] in args['outputs']:
outputname = list(output.keys())[0]
modelresp = np.array(output[outputname])
# Calculate sum of differences
with np.errstate(divide='ignore'): # Ignore warning from taking a log of any zero values
tmp = 20 * np.log10(np.abs(modelresp - refresp) / np.amax(np.abs(refresp)))
# Replace any NaNs or Infs from zero division
tmp[np.invert(np.isfinite(tmp))] = 0
tmp = np.abs(np.sum(tmp)) / len(tmp)
diffdB += tmp
outputs += 1
# Check in case no outputs where found
if outputs == 0:
raise GeneralError('No outputs matching {} were found'.format(args['outputs']))
return diffdB / outputs
def compactness(filename, args):
"""A measure of the compactness of a time domain signal.
Args:
filename (str): Name of output file
args (dict): 'outputs' key with a list of names (IDs) of outputs (rxs) from input file
Returns:
compactness (float): Compactness value of signal
"""
f = h5py.File(filename, 'r')
nrx = f.attrs['nrx']
dt = f.attrs['dt']
iterations = f.attrs['Iterations']
time = np.linspace(0, 1, iterations)
time *= (iterations * dt)
outputsused = False
for rx in range(1, nrx + 1):
output = f['/rxs/rx' + str(rx) + '/']
if output.attrs['Name'] in args['outputs']:
outputname = list(output.keys())[0]
outputdata = output[outputname]
# Get absolute maximum value in signal
peak = np.amax([np.amax(outputdata), np.abs(np.amin(outputdata))])
# Get peaks and troughs in signal
delta = peak / 50 # Considered a peak/trough if it has the max/min value, and was preceded (to the left) by a value lower by delta
maxtab, mintab = peakdet(outputdata, delta)
peaks = maxtab + mintab
peaks.sort()
# Remove any peaks smaller than a threshold
thresholdpeak = 1e-3
peaks = [peak for peak in peaks if np.abs(outputdata[peak]) > thresholdpeak]
# Amplitude ratio of the 1st to 3rd peak - hopefully be a measure of a compact envelope
compactness = np.abs(outputdata[peaks[0]]) / np.abs(outputdata[peaks[2]])
# Percentage of maximum value to measure compactness of signal
# durationthreshold = 2
# Check if there is a peak/trough smaller than threshold
# durationthresholdexist = np.where(np.abs(outputdata[peaks]) < (peak * (durationthreshold / 100)))[0]
# if durationthresholdexist.size == 0:
# compactness = time[peaks[-1]]
# else:
# time2threshold = time[peaks[durationthresholdexist[0]]]
# compactness = time2threshold - time[min(peaks)]
# Check in case no outputs where found
if not outputsused:
raise GeneralError('No outputs matching {} were found'.format(args['outputs']))
return compactness
######################################
# Helper methods for signal analyses #
######################################
def spectral_centroid(x, samplerate):
"""Calculate the spectral centroid of a signal.
Args:
x (float): 1D array containing time domain signal
samplerate (float): Sample rate of signal (Hz)
Returns:
centroid (float): Weighted mean of the frequencies present in the signal
"""
magnitudes = np.abs(np.fft.rfft(x))
length = len(x)
# Positive frequencies
freqs = np.abs(np.fft.fftfreq(length, 1.0 / samplerate)[:length // 2 + 1])
centroid = np.sum(magnitudes * freqs) / np.sum(magnitudes)
return centroid
def zero_crossings(x):
"""Find location of zero crossings in 1D data array.
Args:
x (float): 1D array
Returns:
indexzeros (int): Array of indices of zero crossing locations
"""
pos = x > 0
npos = ~pos
indexzeros = ((pos[:-1] & npos[1:]) | (npos[:-1] & pos[1:])).nonzero()[0]
return indexzeros
def peakdet(v, delta, x=None):
"""Detect peaks and troughs in a vector (adapted from MATLAB script at http://billauer.co.il/peakdet.html).
A point is considered a maximum peak if it has the maximal value, and was preceded (to the left) by a value lower by delta.
Eli Billauer, 3.4.05 (Explicitly not copyrighted).
This function is released to the public domain; Any use is allowed.
Args:
v (float): 1D array
delta (float): threshold for determining peaks/troughs
Returns:
maxtab, mintab (list): Indices of peak/trough locations
"""
maxtab = []
mintab = []
if x is None:
x = np.arange(len(v))
v = np.asarray(v)
if len(v) != len(x):
raise GeneralError('Input vectors v and x must have same length')
if not np.isscalar(delta):
raise GeneralError('Input argument delta must be a scalar')
if delta <= 0:
raise GeneralError('Input argument delta must be positive')
mn, mx = np.Inf, -np.Inf
mnpos, mxpos = np.NaN, np.NaN
lookformax = True
for i in np.arange(len(v)):
this = v[i]
if this > mx:
mx = this
mxpos = x[i]
if this < mn:
mn = this
mnpos = x[i]
if lookformax:
if this < mx - delta:
if int(mxpos) != 0:
maxtab.append(int(mxpos))
mn = this
mnpos = x[i]
lookformax = False
else:
if this > mn + delta:
if int(mnpos) != 0:
mintab.append(int(mnpos))
mx = this
mxpos = x[i]
lookformax = True
return maxtab, mintab

38
user_libs/optimisation_taguchi/plot_results.py
查看文件

@@ -1,38 +0,0 @@
# Copyright (C) 2015-2019, Craig Warren
#
# This module is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License.
# To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/.
#
# Please use the attribution at http://dx.doi.org/10.1190/1.3548506
import argparse
import os
import pickle
from gprMax.optimisation_taguchi import plot_optimisation_history
"""Plots the results (pickled to file) from a Taguchi optimisation process."""
# Parse command line arguments
parser = argparse.ArgumentParser(description='Plots the results (pickled to file) from a Taguchi optimisation process.', usage='cd gprMax; python -m user_libs.optimisation_taguchi.plot_results picklefile')
parser.add_argument('picklefile', help='name of file including path')
args = parser.parse_args()
f = open(args.picklefile, 'rb')
optparamshist = pickle.load(f)
fitnessvalueshist = pickle.load(f)
optparamsinit = pickle.load(f)
print('Optimisations summary for: {}'.format(os.path.split(args.picklefile)[1]))
print('Number of iterations: {:g}'.format(len(fitnessvalueshist)))
print('History of fitness values: {}'.format(fitnessvalueshist))
print('Initial parameter values:')
for item in optparamsinit:
print('\t{}: {}'.format(item[0], item[1]))
print('History of parameter values:')
for key, value in optparamshist.items():
print('\t{}: {}'.format(key, value))
# Plot the history of fitness values and each optimised parameter values for the optimisation
plot_optimisation_history(fitnessvalueshist, optparamshist, optparamsinit)