fixed merge

gprMax/_version.py

@@ -1,4 +1,4 @@
 # This is where the version number is set and read by setup.py and conf.py (for the docs)
 
-__version__ = '3.0.0b26'
+__version__ = '3.0.0b29'
 

gprMax/fractals.py

@@ -23,7 +23,7 @@ from gprMax.constants import floattype, complextype
 from gprMax.utilities import round_value
 
 
-class FractalSurface:
+class FractalSurface(object):
     """Fractal surfaces."""
     
     surfaceIDs = ['xminus', 'xplus', 'yminus', 'yplus', 'zminus', 'zplus']
@@ -106,7 +106,7 @@ class FractalSurface:
         self.fractalsurface = self.fractalsurface * ((self.fractalrange[1] - self.fractalrange[0])/fractalrange) + self.fractalrange[0] - ((self.fractalrange[1] - self.fractalrange[0])/fractalrange) * fractalmin
 
 
-class FractalVolume:
+class FractalVolume(object):
     """Fractal volumes."""
     
     def __init__(self, xs, xf, ys, yf, zs, zf, dimension):
@@ -189,7 +189,7 @@ class FractalVolume:
         self.mask[maskxs:maskxf, maskys:maskyf, maskzs:maskzf] = 1
 
 
-class Grass:
+class Grass(object):
     """Geometry information for blades of grass."""
     
     def __init__(self, numblades):

gprMax/geometry_primitives.pyx

@@ -668,15 +668,33 @@ cpdef void build_voxels_from_array(int xs, int ys, int zs, np.uint16_t[:, :, ::1
     """
 
     cdef Py_ssize_t i, j, k
-    cdef int nx, ny, nz, numID
+    cdef int xf, yf, zf, numID
     
-    nx = data.shape[0]
-    ny = data.shape[1]
-    nz = data.shape[2]
+    # Set bounds to domain if they outside
+    if xs < 0:
+        xs = 0
+    if xs + data.shape[0] >= solid.shape[0]:
+        xf = solid.shape[0] - 1
+    else:
+        xf = xs + data.shape[0]
 
-    for i in range(nx):
-        for j in range(ny):
-            for k in range(nz):
-                numID = data[i, j, k]
-                build_voxel(i + xs, j + ys, k + zs, numID, numID, numID, numID, False, solid, rigidE, rigidH, ID)
+    if ys < 0:
+        ys = 0
+    if ys + data.shape[1] >= solid.shape[1]:
+        yf = solid.shape[1] - 1
+    else:
+        yf = ys + data.shape[1]
+
+    if zs < 0:
+        zs = 0
+    if zs + data.shape[2] >= solid.shape[2]:
+        zf = solid.shape[2] - 1
+    else:
+        zf = zs + data.shape[2]
+
+    for i in range(xs, xf):
+        for j in range(ys, yf):
+            for k in range(zs, zf):
+                numID = data[i - xs, j - ys, k - zs]
+                build_voxel(i, j, k, numID, numID, numID, numID, False, solid, rigidE, rigidH, ID)
 

gprMax/geometry_views.py

@@ -16,7 +16,7 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
-import sys
+import os, sys
 import numpy as np
 from struct import pack
 from gprMax.xdmf import write_output_file
@@ -24,7 +24,7 @@ from gprMax.xdmf import write_output_file
 from gprMax.utilities import round_value
 
 
-class GeometryView:
+class GeometryView(object):
     """Views of the geometry of the model."""
 
     if sys.byteorder == 'little':
@@ -53,7 +53,7 @@ class GeometryView:
         self.dx = dx
         self.dy = dy
         self.dz = dz
-        self.filename = filename
+        self.basefilename = filename
         self.type = type
 
     def write_xdmf(self, modelrun, numbermodelruns, G):
@@ -73,9 +73,9 @@ class GeometryView:
 
         # Construct filename from user-supplied name and model run number
         if numbermodelruns == 1:
-            self.filename = G.inputdirectory + self.filename
+            self.filename = os.path.abspath(os.path.join(G.inputdirectory, self.basefilename))
         else:
-            self.filename = G.inputdirectory + self.filename + str(modelrun)
+            self.filename = os.path.abspath(os.path.join(G.inputdirectory, self.basefilename + str(modelrun)))
 
         if self.type == 'n':
             self.filename += '.vti'
@@ -91,14 +91,18 @@ class GeometryView:
             self.vtk_nycells = self.vtk_yfcells - self.vtk_yscells
             self.vtk_nzcells = self.vtk_zfcells - self.vtk_zscells
 
-            # Create an array and add numeric IDs for PML, sources and receivers
-            self.srcs_rxs_pml = np.zeros((G.nx + 1, G.ny + 1, G.nz + 1), dtype=np.int8)
+            # Create arrays and add numeric IDs for PML, sources and receivers (0 is not set, 1 is PML, srcs and rxs numbered thereafter)
+            self.srcs_pml = np.zeros((G.nx + 1, G.ny + 1, G.nz + 1), dtype=np.int8)
+            self.rxs = np.zeros((G.nx + 1, G.ny + 1, G.nz + 1), dtype=np.int8)
             for pml in G.pmls:
-                self.srcs_rxs_pml[pml.xs:pml.xf, pml.ys:pml.yf, pml.zs:pml.zf] = 1
-            for index, srcrx in enumerate(G.rxs + G.hertziandipoles + G.magneticdipoles + G.voltagesources + G.transmissionlines):
-                self.srcs_rxs_pml[srcrx.xcoord, srcrx.ycoord, srcrx.zcoord] = index + 2
+                self.srcs_pml[pml.xs:pml.xf, pml.ys:pml.yf, pml.zs:pml.zf] = 1
+            for index, src in enumerate(G.hertziandipoles + G.magneticdipoles + G.voltagesources + G.transmissionlines):
+                self.srcs_pml[src.xcoord, src.ycoord, src.zcoord] = index + 2
+            for index, rx in enumerate(G.rxs):
+                self.rxs[rx.xcoord, rx.ycoord, rx.zcoord] = index + 1
 
-            vtk_srcs_rxs_pml_offset = round_value((np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + np.dtype(np.uint32).itemsize)
+            vtk_srcs_pml_offset = round_value((np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + np.dtype(np.uint32).itemsize)
+            vtk_rxs_offset = round_value((np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + np.dtype(np.uint32).itemsize + (np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + np.dtype(np.uint32).itemsize)
 
             with open(self.filename, 'wb') as f:
                 f.write('<?xml version="1.0"?>\n'.encode('utf-8'))
@@ -107,7 +111,8 @@ class GeometryView:
                 f.write('<Piece Extent="{} {} {} {} {} {}">\n'.format(self.vtk_xscells, self.vtk_xfcells, self.vtk_yscells, self.vtk_yfcells, self.vtk_zscells, self.vtk_zfcells).encode('utf-8'))
                 f.write('<CellData Scalars="Material">\n'.encode('utf-8'))
                 f.write('<DataArray type="UInt32" Name="Material" format="appended" offset="0" />\n'.encode('utf-8'))
-                f.write('<DataArray type="Int8" Name="Sources_Receivers_PML" format="appended" offset="{}" />\n'.format(vtk_srcs_rxs_pml_offset).encode('utf-8'))
+                f.write('<DataArray type="Int8" Name="Sources_PML" format="appended" offset="{}" />\n'.format(vtk_srcs_pml_offset).encode('utf-8'))
+                f.write('<DataArray type="Int8" Name="Receivers" format="appended" offset="{}" />\n'.format(vtk_rxs_offset).encode('utf-8'))
                 f.write('</CellData>\n'.encode('utf-8'))
                 f.write('</Piece>\n</ImageData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
 
@@ -120,13 +125,21 @@ class GeometryView:
                         for i in range(self.xs, self.xf, self.dx):
                             f.write(pack('I', G.solid[i, j, k]))
 
-                # Write source/receiver IDs
+                # Write source/PML IDs
                 datasize = int(np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)
                 f.write(pack('I', datasize))
                 for k in range(self.zs, self.zf, self.dz):
                     for j in range(self.ys, self.yf, self.dy):
                         for i in range(self.xs, self.xf, self.dx):
-                            f.write(pack('b', self.srcs_rxs_pml[i, j, k]))
+                            f.write(pack('b', self.srcs_pml[i, j, k]))
+
+                # Write receiver IDs
+                datasize = int(np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells)
+                f.write(pack('I', datasize))
+                for k in range(self.zs, self.zf, self.dz):
+                    for j in range(self.ys, self.yf, self.dy):
+                        for i in range(self.xs, self.xf, self.dx):
+                            f.write(pack('b', self.rxs[i, j, k]))
 
                 f.write('\n</AppendedData>\n</VTKFile>'.encode('utf-8'))
 
@@ -244,8 +257,8 @@ class GeometryView:
             f.write('<Material name="{}">{}</Material>\n'.format(material.ID, material.numID).encode('utf-8'))
         if not materialsonly:
             f.write('<PML name="PML boundary region">1</PML>\n'.encode('utf-8'))
-            for index, srcrx in enumerate(G.rxs + G.hertziandipoles + G.magneticdipoles + G.voltagesources + G.transmissionlines):
-                f.write('<Sources_Receivers name="{}">{}</Sources_Receivers>\n'.format(srcrx.ID, index + 2).encode('utf-8'))
+            for index, src in enumerate(G.hertziandipoles + G.magneticdipoles + G.voltagesources + G.transmissionlines):
+                f.write('<Sources name="{}">{}</Sources>\n'.format(src.ID, index + 2).encode('utf-8'))
+            for index, rx in enumerate(G.rxs):
+                f.write('<Receivers name="{}">{}</Receivers>\n'.format(rx.ID, index + 1).encode('utf-8'))
         f.write('</gprMax>\n'.encode('utf-8'))
-
-

gprMax/gprMax.py

@@ -53,15 +53,16 @@ def main():
     parser.add_argument('-mpi', action='store_true', default=False, help='switch on MPI task farm')
     parser.add_argument('-benchmark', action='store_true', default=False, help='switch on benchmarking mode')
     parser.add_argument('--geometry-only', action='store_true', default=False, help='only build model and produce geometry file(s)')
+    parser.add_argument('--geometry-fixed', action='store_true', default=False, help='do not reprocess model geometry for multiple model runs')
     parser.add_argument('--write-processed', action='store_true', default=False, help='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='optimise parameters using the Taguchi optimisation method')
     args = parser.parse_args()
     numbermodelruns = args.n
-    inputdirectory = os.path.dirname(os.path.abspath(args.inputfile)) + os.sep
-    inputfile = inputdirectory + os.path.basename(args.inputfile)
+    inputdirectory = os.path.dirname(os.path.abspath(args.inputfile))
+    inputfile = os.path.abspath(os.path.join(inputdirectory, os.path.basename(args.inputfile)))
 
     # 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': numbermodelruns, 'inputdirectory': inputdirectory}
+    usernamespace = {'c': c, 'e0': e0, 'm0': m0, 'z0': z0, 'number_model_runs': numbermodelruns, 'input_directory': inputdirectory}
 
     # Process for Taguchi optimisation
     if args.opt_taguchi:
@@ -239,107 +240,143 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
     # Monitor memory usage
     p = psutil.Process()
 
-    print('\n{}\n\nModel input file: {}\n'.format(68*'*', inputfile))
+    # Declare variable to hold FDTDGrid class
+    global G
 
-    # Add the current model run to namespace that can be accessed by user in any Python code blocks in input file
-    usernamespace['current_model_run'] = modelrun
-    print('Constants/variables available for Python scripting: {}\n'.format(usernamespace))
+    # Normal model reading/building process; bypassed if geometry information to be reused
+    if not 'G' in globals():
+        print('\n{}\n\nModel input file: {}\n'.format(68*'*', inputfile))
 
-    # Process any user input Python commands
-    processedlines = process_python_include_code(inputfile, usernamespace)
+        # Add the current model run to namespace that can be accessed by user in any Python code blocks in input file
+        usernamespace['current_model_run'] = modelrun
+        print('Constants/variables available for Python scripting: {}\n'.format(usernamespace))
 
-    # Write a file containing the input commands after Python blocks have been processed
-    if args.write_processed:
-        write_processed_file(inputfile, modelrun, numbermodelruns, processedlines)
+        # Read input file and process any Python or include commands
+        processedlines = process_python_include_code(inputfile, usernamespace)
 
-    # Check validity of command names & that essential commands are present
-    singlecmds, multicmds, geometry = check_cmd_names(processedlines)
+        # Write a file containing the input commands after Python or include commands have been processed
+        if args.write_processed:
+            write_processed_file(inputfile, modelrun, numbermodelruns, processedlines)
 
-    # Initialise an instance of the FDTDGrid class
-    G = FDTDGrid()
-    G.inputdirectory = usernamespace['inputdirectory']
+        # Check validity of command names and that essential commands are present
+        singlecmds, multicmds, geometry = check_cmd_names(processedlines)
 
-    # Create built-in materials
-    m = Material(0, 'pec', G)
-    m.average = False
-    G.materials.append(m)
-    m = Material(1, 'free_space', G)
-    G.materials.append(m)
+        # Initialise an instance of the FDTDGrid class
+        G = FDTDGrid()
+        G.inputfilename = os.path.split(inputfile)[1]
+        G.inputdirectory = os.path.dirname(os.path.abspath(inputfile))
 
-    # Process parameters for commands that can only occur once in the model
-    process_singlecmds(singlecmds, G)
+        # Create built-in materials
+        m = Material(0, 'pec', G)
+        m.average = False
+        G.materials.append(m)
+        m = Material(1, 'free_space', G)
+        G.materials.append(m)
 
-    # Process parameters for commands that can occur multiple times in the model
-    process_multicmds(multicmds, G)
+        # Process parameters for commands that can only occur once in the model
+        process_singlecmds(singlecmds, G)
 
-    # Initialise an array for volumetric material IDs (solid), boolean arrays for specifying materials not to be averaged (rigid),
-    # an array for cell edge IDs (ID), and arrays for the field components.
-    G.initialise_std_arrays()
+        # Process parameters for commands that can occur multiple times in the model
+        process_multicmds(multicmds, G)
 
-    # Process the geometry commands in the order they were given
-    tinputprocstart = perf_counter()
-    process_geometrycmds(geometry, G)
-    tinputprocend = perf_counter()
-    print('\nInput file processed in [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=int(tinputprocend - tinputprocstart))))
+        # Initialise an array for volumetric material IDs (solid), boolean arrays for specifying materials not to be averaged (rigid),
+        # an array for cell edge IDs (ID)
+        G.initialise_geometry_arrays()
 
-    # Build the PML and calculate initial coefficients
-    build_pmls(G)
+        # Initialise arrays for the field components
+        G.initialise_field_arrays()
 
-    # Build the model, i.e. set the material properties (ID) for every edge of every Yee cell
-    tbuildstart = perf_counter()
-    build_electric_components(G.solid, G.rigidE, G.ID, G)
-    build_magnetic_components(G.solid, G.rigidH, G.ID, G)
-    tbuildend = perf_counter()
-    print('\nModel built in [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=int(tbuildend - tbuildstart))))
+        # Process geometry commands in the order they were given
+        tinputprocstart = perf_counter()
+        process_geometrycmds(geometry, G)
+        tinputprocend = perf_counter()
+        print('\nInput file processed in [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=int(tinputprocend - tinputprocstart))))
 
-    # Process any voltage sources (that have resistance) to create a new material at the source location
-    for voltagesource in G.voltagesources:
-        voltagesource.create_material(G)
+        # Build the PML and calculate initial coefficients
+        build_pmls(G)
 
-    # Initialise arrays of update coefficients to pass to update functions
-    G.initialise_std_updatecoeff_arrays()
+        # Build the model, i.e. set the material properties (ID) for every edge of every Yee cell
+        tbuildstart = perf_counter()
+        build_electric_components(G.solid, G.rigidE, G.ID, G)
+        build_magnetic_components(G.solid, G.rigidH, G.ID, G)
+        tbuildend = perf_counter()
+        print('\nModel built in [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=int(tbuildend - tbuildstart))))
 
-    # Initialise arrays of update coefficients and temporary values if there are any dispersive materials
-    if Material.maxpoles != 0:
-        G.initialise_dispersive_arrays()
+        # Process any voltage sources (that have resistance) to create a new material at the source location
+        for voltagesource in G.voltagesources:
+            voltagesource.create_material(G)
 
-    # Calculate update coefficients, store in arrays, and list materials in model
-    if G.messages:
-        print('\nMaterials:\n')
-        print('ID\tName\t\tProperties')
-        print('{}'.format('-'*50))
-    for material in G.materials:
+        # Initialise arrays of update coefficients to pass to update functions
+        G.initialise_std_update_coeff_arrays()
 
-        # Calculate update coefficients for material
-        material.calculate_update_coeffsE(G)
-        material.calculate_update_coeffsH(G)
-
-        # Store all update coefficients together
-        G.updatecoeffsE[material.numID, :] = material.CA, material.CBx, material.CBy, material.CBz, material.srce
-        G.updatecoeffsH[material.numID, :] = material.DA, material.DBx, material.DBy, material.DBz, material.srcm
-
-        # Store coefficients for any dispersive materials
+        # Initialise arrays of update coefficients and temporary values if there are any dispersive materials
         if Material.maxpoles != 0:
-            z = 0
-            for pole in range(Material.maxpoles):
-                G.updatecoeffsdispersive[material.numID, z:z+3] = e0 * material.eqt2[pole], material.eqt[pole], material.zt[pole]
-                z += 3
+            G.initialise_dispersive_arrays()
 
+        # Calculate update coefficients, store in arrays, and list materials in model
         if G.messages:
-            if material.deltaer and material.tau:
-                tmp = 'delta_epsr={}, tau={} secs; '.format(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer), ', '.join('{:g}'.format(tau) for tau in material.tau))
-            else:
-                tmp = ''
-            if material.average:
-                dielectricsmoothing = 'dielectric smoothing permitted.'
-            else:
-                dielectricsmoothing = 'dielectric smoothing not permitted.'
-            print('{:3}\t{:12}\tepsr={:g}, sig={:g} S/m; mur={:g}, sig*={:g} S/m; '.format(material.numID, material.ID, material.er, material.se, material.mr, material.sm) + tmp + dielectricsmoothing)
+            print('\nMaterials:\n')
+            print('ID\tName\t\tProperties')
+            print('{}'.format('-'*50))
+        for material in G.materials:
 
-    # Check to see if numerical dispersion might be a problem
-    resolution = dispersion_check(G)
-    if resolution != 0 and max((G.dx, G.dy, G.dz)) > resolution:
-        print('\nWARNING: Potential numerical dispersion in the simulation. Check the spatial discretisation against the smallest wavelength present. Suggested resolution should be less than {:g}m'.format(resolution))
+            # Calculate update coefficients for material
+            material.calculate_update_coeffsE(G)
+            material.calculate_update_coeffsH(G)
+
+            # Store all update coefficients together
+            G.updatecoeffsE[material.numID, :] = material.CA, material.CBx, material.CBy, material.CBz, material.srce
+            G.updatecoeffsH[material.numID, :] = material.DA, material.DBx, material.DBy, material.DBz, material.srcm
+
+            # Store coefficients for any dispersive materials
+            if Material.maxpoles != 0:
+                z = 0
+                for pole in range(Material.maxpoles):
+                    G.updatecoeffsdispersive[material.numID, z:z+3] = e0 * material.eqt2[pole], material.eqt[pole], material.zt[pole]
+                    z += 3
+
+            if G.messages:
+                if material.deltaer and material.tau:
+                    tmp = 'delta_epsr={}, tau={} secs; '.format(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer), ', '.join('{:g}'.format(tau) for tau in material.tau))
+                else:
+                    tmp = ''
+                if material.average:
+                    dielectricsmoothing = 'dielectric smoothing permitted.'
+                else:
+                    dielectricsmoothing = 'dielectric smoothing not permitted.'
+                print('{:3}\t{:12}\tepsr={:g}, sig={:g} S/m; mur={:g}, sig*={:g} S/m; '.format(material.numID, material.ID, material.er, material.se, material.mr, material.sm) + tmp + dielectricsmoothing)
+
+        # Check to see if numerical dispersion might be a problem
+        resolution = dispersion_check(G)
+        if resolution and max((G.dx, G.dy, G.dz)) > resolution:
+            print('\nWARNING: Potential numerical dispersion in the simulation. Check the spatial discretisation against the smallest wavelength present. Suggested resolution should be less than {:g}m'.format(resolution))
+
+    # If geometry information to be reused between model runs
+    else:
+        # Clear arrays for field components
+        G.initialise_field_arrays()
+
+        # Clear arrays for fields in PML
+        for pml in G.pmls:
+            pml.initialise_field_arrays()
+
+    # Adjust position of simple sources and receivers if required
+    if G.srcstepx > 0 or G.srcstepy > 0 or G.srcstepz > 0:
+        for source in itertools.chain(G.hertziandipoles, G.magneticdipoles):
+            if modelrun == 1:
+                if source.xcoord + G.srcstepx * (numbermodelruns - 1) > G.nx or source.ycoord + G.srcstepy * (numbermodelruns - 1) > G.ny or source.zcoord + G.srcstepz * (numbermodelruns - 1) > G.nz:
+                    raise GeneralError('Source(s) will be stepped to a position outside the domain.')
+            source.xcoord = source.xcoordbase + (modelrun - 1) * G.srcstepx
+            source.ycoord = source.ycoordbase + (modelrun - 1) * G.srcstepy
+            source.zcoord = source.zcoordbase + (modelrun - 1) * G.srcstepz
+    if G.rxstepx > 0 or G.rxstepy > 0 or G.rxstepz > 0:
+        for receiver in G.rxs:
+            if modelrun == 1:
+                if receiver.xcoord + G.rxstepx * (numbermodelruns - 1) > G.nx or receiver.ycoord + G.rxstepy * (numbermodelruns - 1) > G.ny or receiver.zcoord + G.rxstepz * (numbermodelruns - 1) > G.nz:
+                    raise GeneralError('Receiver(s) will be stepped to a position outside the domain.')
+            receiver.xcoord = receiver.xcoordbase + (modelrun - 1) * G.rxstepx
+            receiver.ycoord = receiver.ycoordbase + (modelrun - 1) * G.rxstepy
+            receiver.zcoord = receiver.zcoordbase + (modelrun - 1) * G.rxstepz
 
     # Write files for any geometry views
     if not G.geometryviews and args.geometry_only:
@@ -352,31 +389,13 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
         tgeoend = perf_counter()
         print('\nGeometry file(s) written in [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=int(tgeoend - tgeostart))))
 
-    # Run simulation if not doing only geometry
+    # Run simulation (if not doing geometry only)
     if not args.geometry_only:
 
         # Prepare any snapshot files
         for snapshot in G.snapshots:
             snapshot.prepare_vtk_imagedata(modelrun, numbermodelruns, G)
 
-        # Adjust position of sources and receivers if required
-        if G.srcstepx > 0 or G.srcstepy > 0 or G.srcstepz > 0:
-            for source in itertools.chain(G.hertziandipoles, G.magneticdipoles, G.voltagesources, G.transmissionlines):
-                if modelrun == 1:
-                    if source.xcoord + G.srcstepx * (numbermodelruns - 1) > G.nx or source.ycoord + G.srcstepy * (numbermodelruns - 1) > G.ny or source.zcoord + G.srcstepz * (numbermodelruns - 1) > G.nz:
-                        raise GeneralError('Source(s) will be stepped to a position outside the domain.')
-                source.xcoord += (modelrun - 1) * G.srcstepx
-                source.ycoord += (modelrun - 1) * G.srcstepy
-                source.zcoord += (modelrun - 1) * G.srcstepz
-        if G.rxstepx > 0 or G.rxstepy > 0 or G.rxstepz > 0:
-            for receiver in G.rxs:
-                if modelrun == 1:
-                    if receiver.xcoord + G.rxstepx * (numbermodelruns - 1) > G.nx or receiver.ycoord + G.rxstepy * (numbermodelruns - 1) > G.ny or receiver.zcoord + G.rxstepz * (numbermodelruns - 1) > G.nz:
-                        raise GeneralError('Receiver(s) will be stepped to a position outside the domain.')
-                receiver.xcoord += (modelrun - 1) * G.rxstepx
-                receiver.ycoord += (modelrun - 1) * G.rxstepy
-                receiver.zcoord += (modelrun - 1) * G.rxstepz
-
         # Prepare output file
         inputfileparts = os.path.splitext(inputfile)
         if numbermodelruns == 1:
@@ -466,6 +485,10 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
         print('\n\nSolving took [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=int(tsolveend - tsolvestart))))
         print('Peak memory (approx) used: {}'.format(human_size(p.memory_info().rss)))
 
+        # If geometry information to be reused between model runs then FDTDGrid class instance must be global so that it persists
+        if not args.geometry_fixed:
+            del G
+
         return int(tsolveend - tsolvestart)
 
 

gprMax/grid.py

@@ -22,6 +22,7 @@ import matplotlib.pyplot as plt
 from gprMax.constants import c, floattype, complextype
 from gprMax.materials import Material
 
+
 class Grid():
 
     def __init__(self, grid):
@@ -46,10 +47,12 @@ class Grid():
     def get(self, i, j, k):
         return self.grid[i, j, k]
 
+
 class FDTDGrid(Grid):
     """Holds attributes associated with the entire grid. A convenient way for accessing regularly used parameters."""
 
     def __init__(self):
+        self.inputfilename = ''
         self.inputdirectory = ''
         self.title = ''
         self.messages = True
@@ -76,25 +79,27 @@ class FDTDGrid(Grid):
         self.hertziandipoles = []
         self.magneticdipoles = []
         self.transmissionlines = []
+        self.rxs = []
         self.srcstepx = 0
         self.srcstepy = 0
         self.srcstepz = 0
         self.rxstepx = 0
         self.rxstepy = 0
         self.rxstepz = 0
-        self.rxs = []
         self.snapshots = []
 
-    def initialise_std_arrays(self):
+    def initialise_geometry_arrays(self):
         """Initialise an array for volumetric material IDs (solid); boolean arrays for specifying whether materials can have dielectric smoothing (rigid);
-            an array for cell edge IDs (ID); and arrays for the electric and magnetic field components. Solid and ID arrays are initialised to free_space (one); rigid arrays
-            to allow dielectric smoothing (zero).
+            and an array for cell edge IDs (ID). Solid and ID arrays are initialised to free_space (one); rigid arrays to allow dielectric smoothing (zero).
         """
         self.solid = np.ones((self.nx + 1, self.ny + 1, self.nz + 1), dtype=np.uint32)
         self.rigidE = np.zeros((12, self.nx + 1, self.ny + 1, self.nz + 1), dtype=np.int8)
         self.rigidH = np.zeros((6, self.nx + 1, self.ny + 1, self.nz + 1), dtype=np.int8)
         self.IDlookup = {'Ex': 0, 'Ey': 1, 'Ez': 2, 'Hx': 3, 'Hy': 4, 'Hz': 5}
         self.ID = np.ones((6, self.nx + 1, self.ny + 1, self.nz + 1), dtype=np.uint32)
+
+    def initialise_field_arrays(self):
+        """Initialise arrays for the electric and magnetic field components."""
         self.Ex = np.zeros((self.nx, self.ny + 1, self.nz + 1), dtype=floattype)
         self.Ey = np.zeros((self.nx + 1, self.ny, self.nz + 1), dtype=floattype)
         self.Ez = np.zeros((self.nx + 1, self.ny + 1, self.nz), dtype=floattype)
@@ -102,7 +107,7 @@ class FDTDGrid(Grid):
         self.Hy = np.zeros((self.nx, self.ny + 1, self.nz), dtype=floattype)
         self.Hz = np.zeros((self.nx, self.ny, self.nz + 1), dtype=floattype)
 
-    def initialise_std_updatecoeff_arrays(self):
+    def initialise_std_update_coeff_arrays(self):
         """Initialise arrays for storing update coefficients."""
         self.updatecoeffsE = np.zeros((len(self.materials), 5), dtype=floattype)
         self.updatecoeffsH = np.zeros((len(self.materials), 5), dtype=floattype)
@@ -132,31 +137,38 @@ def dispersion_check(G):
     maxfreqs = []
     for waveform in G.waveforms:
 
-        # User-defined waveform
-        if waveform.uservalues is not None:
-            waveformvalues = waveform.uservalues
+        if waveform.type == 'sine' or waveform.type == 'contsine':
+            maxfreqs.append(4 * waveform.freq)
+
+        elif waveform.type =='impulse':
+            pass
 
-        # Built-in waveform
         else:
-            time = np.linspace(0, 1, G.iterations)
-            time *= (G.iterations * G.dt)
-            waveformvalues = np.zeros(len(time))
-            timeiter = np.nditer(time, flags=['c_index'])
+            # User-defined waveform
+            if waveform.type == 'user':
+                waveformvalues = waveform.uservalues
 
-            while not timeiter.finished:
-                waveformvalues[timeiter.index] = waveform.calculate_value(timeiter[0], G.dt)
-                timeiter.iternext()
+            # Built-in waveform
+            else:
+                time = np.linspace(0, 1, G.iterations)
+                time *= (G.iterations * G.dt)
+                waveformvalues = np.zeros(len(time))
+                timeiter = np.nditer(time, flags=['c_index'])
 
-        # Calculate magnitude of frequency spectra of waveform
-        power = 20 * np.log10(np.abs(np.fft.fft(waveformvalues))**2)
-        freqs = np.fft.fftfreq(power.size, d=G.dt)
+                while not timeiter.finished:
+                    waveformvalues[timeiter.index] = waveform.calculate_value(timeiter[0], G.dt)
+                    timeiter.iternext()
 
-        # Shift powers so that frequency with maximum power is at zero decibels
-        power -= np.amax(power)
+            # Calculate magnitude of frequency spectra of waveform
+            power = 20 * np.log10(np.abs(np.fft.fft(waveformvalues))**2)
+            freqs = np.fft.fftfreq(power.size, d=G.dt)
 
-        # Set maximum frequency to -60dB from maximum power
-        freq = np.where((np.amax(power[1::]) - power[1::]) > 60)[0][0] + 1
-        maxfreqs.append(freqs[freq])
+            # Shift powers so that frequency with maximum power is at zero decibels
+            power -= np.amax(power)
+
+            # Set maximum frequency to -60dB from maximum power
+            freq = np.where((np.amax(power[1::]) - power[1::]) > 60)[0][0] + 1
+            maxfreqs.append(freqs[freq])
 
     if maxfreqs:
         maxfreq = max(maxfreqs)
@@ -175,7 +187,7 @@ def dispersion_check(G):
         resolution = minwavelength / resolvedsteps
 
     else:
-        resolution = 0
+        resolution = False
 
     return resolution
 

gprMax/input_cmds_file.py

@@ -32,7 +32,9 @@ def process_python_include_code(inputfile, usernamespace):
     Returns:
         processedlines (list): Input commands after Python processing.
     """
-        
+    
+    userpython = False
+    
     with open(inputfile, 'r') as f:
         # Strip out any newline characters and comments that must begin with double hashes
         inputlines = [line.rstrip() for line in f if(not line.startswith('##') and line.rstrip('\n'))]
@@ -45,6 +47,9 @@ def process_python_include_code(inputfile, usernamespace):
         
         # Process any Python code
         if(inputlines[x].startswith('#python:')):
+            # Save stdout location to restore later
+            stdout = sys.stdout
+            
             # String to hold Python code to be executed
             pythoncode = ''
             x += 1
@@ -67,12 +72,15 @@ def process_python_include_code(inputfile, usernamespace):
             # Add processed Python code to list
             processedlines.extend(codeproc)
     
+            # Reset stdio
+            sys.stdout = stdout
+    
         # Process any include commands
-        elif(inputlines[x].startswith('#include:')):
+        elif(inputlines[x].startswith('#include_file:')):
             includefile = inputlines[x].split()
             
             if len(includefile) != 2:
-                raise CmdInputError('#include requires exactly one parameter')
+                raise CmdInputError('#include_file requires exactly one parameter')
             
             includefile = includefile[1]
             
@@ -94,8 +102,6 @@ def process_python_include_code(inputfile, usernamespace):
             processedlines.append(inputlines[x])
 
         x += 1
-
-    sys.stdout = sys.__stdout__ # Reset stdio
     
     return processedlines
 

gprMax/input_cmds_geometry.py

@@ -53,7 +53,7 @@ def process_geometrycmds(geometry, G):
             
             # See if material file exists at specified path and if not try input file directory
             if not os.path.isfile(matfile):
-                matfile = os.path.join(G.inputdirectory, matfile)
+                matfile = os.path.abspath(os.path.join(G.inputdirectory, matfile))
             
             # Read materials from file
             with open(matfile, 'r') as f:
@@ -70,25 +70,21 @@ def process_geometrycmds(geometry, G):
             
             # See if geometry object file exists at specified path and if not try input file directory
             if not os.path.isfile(geofile):
-                geofile = os.path.join(G.inputdirectory, geofile)
+                geofile = os.path.abspath(os.path.join(G.inputdirectory, geofile))
             
-            # Open geometry object file and read spatial resolution attribute
+            # Open geometry object file and read/check spatial resolution attribute
             f = h5py.File(geofile, 'r')
             dx_dy_dz = f.attrs['dx, dy, dz']
             if dx_dy_dz[0] != G.dx or dx_dy_dz[1] != G.dy or dx_dy_dz[2] != G.dz:
-                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires the spatial resolution of the objects file to match the spatial resolution of the model')
+                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires the spatial resolution of the geometry objects file to match the spatial resolution of the model')
 
             data = f['/data'][:]
-            nx = data.shape[0]
-            ny = data.shape[1]
-            nz = data.shape[2]
-            
-            if (xs + nx) > G.nx or (ys + ny) > G.ny or (zs + nz) > G.nz:
-                raise CmdInputError("'" + ' '.join(tmp) + "'" + ' the requested geometry objects do not fit within the model domain')
-        
             data += numexistmaterials
             build_voxels_from_array(xs, ys, zs, data, G.solid, G.rigidE, G.rigidH, G.ID)
     
+            if G.messages:
+                print('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))
+    
         
         elif tmp[0] == '#edge:':
             if len(tmp) != 8:

gprMax/input_cmds_multiuse.py

@@ -160,6 +160,9 @@ def process_multicmds(multicmds, G):
             h.xcoord = xcoord
             h.ycoord = ycoord
             h.zcoord = zcoord
+            h.xcoordbase = xcoord
+            h.ycoordbase = ycoord
+            h.zcoordbase = zcoord
             h.ID = 'HertzianDipole(' + str(h.xcoord) + ',' + str(h.ycoord) + ',' + str(h.zcoord) + ')'
             h.waveformID = tmp[4]
             
@@ -222,6 +225,9 @@ def process_multicmds(multicmds, G):
             m.xcoord = xcoord
             m.ycoord = ycoord
             m.zcoord = zcoord
+            m.xcoordbase = xcoord
+            m.ycoordbase = ycoord
+            m.zcoordbase = zcoord
             m.ID = 'MagneticDipole(' + str(m.xcoord) + ',' + str(m.ycoord) + ',' + str(m.zcoord) + ')'
             m.waveformID = tmp[4]
             
@@ -340,7 +346,13 @@ def process_multicmds(multicmds, G):
             if xcoord < G.pmlthickness[0] or xcoord > G.nx - G.pmlthickness[3] or ycoord < G.pmlthickness[1] or ycoord > G.ny - G.pmlthickness[4] or zcoord < G.pmlthickness[2] or zcoord > G.nz - G.pmlthickness[5]:
                 print("WARNING: '" + cmdname + ': ' + ' '.join(tmp) + "'" + ' sources and receivers should not normally be positioned within the PML.\n')
             
-            r = Rx(xcoord=xcoord, ycoord=ycoord, zcoord=zcoord)
+            r = Rx()
+            r.xcoord = xcoord
+            r.ycoord = ycoord
+            r.zcoord = zcoord
+            r.xcoordbase = xcoord
+            r.ycoordbase = ycoord
+            r.zcoordbase = zcoord
             
             # If no ID or outputs are specified, use default i.e Ex, Ey, Ez, Hx, Hy, Hz, Ix, Iy, Iz
             if len(tmp) == 3:
@@ -403,7 +415,13 @@ def process_multicmds(multicmds, G):
             for x in range(xs, xf, dx):
                 for y in range(ys, yf, dy):
                     for z in range(zs, zf, dz):
-                        r = Rx(xcoord=x, ycoord=y, zcoord=z)
+                        r = Rx()
+                        r.xcoord = x
+                        r.ycoord = y
+                        r.zcoord = z
+                        r.xcoordbase = x
+                        r.ycoordbase = y
+                        r.zcoordbase = z
                         r.ID = 'Rx(' + str(x) + ',' + str(y) + ',' + str(z) + ')'
                         G.rxs.append(r)
 
@@ -429,15 +447,16 @@ def process_multicmds(multicmds, G):
             dy = round_value(float(tmp[7])/G.dy)
             dz = round_value(float(tmp[8])/G.dz)
             
+            # If number of iterations given
+            try:
+                time = int(tmp[9])
             # If real floating point value given
-            if '.' in tmp[9] or 'e' in tmp[9]:
-                if float(tmp[9]) > 0:
-                    time = round_value((float(tmp[9]) / G.dt)) + 1
+            except:
+                time = float(tmp[9])
+                if time > 0:
+                    time = round_value((time / G.dt)) + 1
                 else:
                     raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' time value must be greater than zero')
-            # If number of iterations given
-            else:
-                time = int(tmp[9])
             
             if xs < 0 or xs > G.nx:
                 raise CmdInputError("'" + cmdname + ': ' + ' '.join(tmp) + "'" + ' the lower x-coordinate {:g}m is not within the model domain'.format(xs * G.dx))
@@ -463,7 +482,7 @@ def process_multicmds(multicmds, G):
             s = Snapshot(xs, ys, zs, xf, yf, zf, dx, dy, dz, time, tmp[10])
 
             if G.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, dx * G.dy, dx * G.dz, s.time * G.dt, s.filename))
+                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, dx * G.dy, dx * G.dz, s.time * G.dt, s.basefilename))
                     
             G.snapshots.append(s)
 
@@ -533,7 +552,7 @@ def process_multicmds(multicmds, G):
                     Material.maxpoles = material.poles
 
                 if G.messages:
-                    print('Debye-type disperion added to {} with delta_epsr={}, 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)))
+                    print('Debye disperion added to {} with delta_epsr={}, 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'
     if multicmds[cmdname] != 'None':
@@ -569,7 +588,7 @@ def process_multicmds(multicmds, G):
                     Material.maxpoles = material.poles
 
                 if G.messages:
-                    print('Lorentz-type disperion added to {} with delta_epsr={}, 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)))
+                    print('Lorentz disperion added to {} with delta_epsr={}, 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'
@@ -605,7 +624,7 @@ def process_multicmds(multicmds, G):
                     Material.maxpoles = material.poles
 
                 if G.messages:
-                    print('Drude-type 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)))
+                    print('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'
@@ -685,7 +704,7 @@ def process_multicmds(multicmds, G):
             g = GeometryView(xs, ys, zs, xf, yf, zf, dx, dy, dz, tmp[9], tmp[10].lower())
 
             if G.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, 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, g.filename))
+                print('Geometry view from {:g}m, {:g}m, {:g}m, to {:g}m, {:g}m, {:g}m, discretisation {:g}m, {:g}m, {:g}m, 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, g.basefilename))
 
             # Append the new GeometryView object to the geometry views list
             G.geometryviews.append(g)

gprMax/input_cmds_singleuse.py

@@ -164,17 +164,20 @@ def process_singlecmds(singlecmds, G):
     if len(tmp) != 1:
         raise CmdInputError(cmd + ' requires exactly one parameter to specify the time window. Either in seconds or number of iterations.')
     tmp = tmp[0].lower()
+
+    # If number of iterations given
+    try:
+        tmp = int(tmp)
+        G.timewindow = (tmp - 1) * G.dt
+        G.iterations = tmp
     # If real floating point value given
-    if '.' in tmp or 'e' in tmp:
-        if float(tmp) > 0:
-            G.timewindow = float(tmp)
-            G.iterations = round_value((float(tmp) / G.dt)) + 1
+    except:
+        tmp = float(tmp)
+        if tmp > 0:
+            G.timewindow = tmp
+            G.iterations = round_value((tmp / G.dt)) + 1
         else:
             raise CmdInputError(cmd + ' must have a value greater than zero')
-    # If number of iterations given
-    else:
-        G.timewindow = (int(tmp) - 1) * G.dt
-        G.iterations = int(tmp)
     if G.messages:
         print('Time window: {:g} secs ({} iterations)'.format(G.timewindow, G.iterations))
 
@@ -203,7 +206,7 @@ def process_singlecmds(singlecmds, G):
         G.srcstepy = round_value(float(tmp[1])/G.dy)
         G.srcstepz = round_value(float(tmp[2])/G.dz)
         if G.messages:
-            print('All sources will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.srcstepx * G.dx, G.srcstepy * G.dy, G.srcstepz * G.dz))
+            print('Simple sources will step {:g}m, {:g}m, {:g}m for each model run.'.format(G.srcstepx * G.dx, G.srcstepy * G.dy, G.srcstepz * G.dz))
 
 
     # rx_steps
@@ -229,7 +232,7 @@ def process_singlecmds(singlecmds, G):
 
         # See if file exists at specified path and if not try input file directory
         if not os.path.isfile(excitationfile):
-            excitationfile = os.path.join(G.inputdirectory, excitationfile)
+            excitationfile = os.path.abspath(os.path.join(G.inputdirectory, excitationfile))
 
         # Get waveform names
         with open(excitationfile, 'r') as f:

gprMax/materials.py

@@ -21,7 +21,7 @@ import numpy as np
 from gprMax.constants import e0, m0, floattype, complextype
 
 
-class Material():
+class Material(object):
     """Materials, their properties and update coefficients."""
     
     # Maximum number of dispersive material poles in a model
@@ -144,7 +144,7 @@ class Material():
             self.srce = 1 / EA
 
 
-class PeplinskiSoil:
+class PeplinskiSoil(object):
     """Soil objects that are characterised according to a mixing model by Peplinski (http://dx.doi.org/10.1109/36.387598)."""
     
     def __init__(self, ID, sandfraction, clayfraction, bulkdensity, sandpartdensity, watervolfraction):

gprMax/optimisation_taguchi.py

@@ -62,11 +62,17 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace):
     optparamshist = OrderedDict((key, list()) for key in optparams)
     
     # Import specified fitness function
-    fitness_metric = getattr(importlib.import_module('user_libs.optimisation_taguchi.optimisation_taguchi_fitness'), fitness['name'])
+    fitness_metric = getattr(importlib.import_module('user_libs.optimisation_taguchi.fitness_functions'), fitness['name'])
 
     # Select OA
     OA, N, cols, k, s, t = construct_OA(optparams)
-    print('\n{}\n\nTaguchi optimisation: orthogonal array with {} experiments, {} parameters ({} used), {} levels, and strength {} will be used.'.format(68*'*', N, cols, k, s, t))
+    print('\n{}\nTaguchi optimisation...\n'.format(68*'*'))
+    print('\tOrthogonal 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('\tParameters to optimise with ranges: {}'.format(str(tmp).strip('[]')))
+    print('\tOutput name(s) from model: {}'.format(fitness['args']['outputs']))
+    print('\tFitness function {} with stopping criterion {:g}'.format(fitness['name'], fitness['stop']))
+    print('\tMaximum iterations: {:g}'.format(maxiterations))
     
     # Initialise arrays and lists to store parameters required throughout optimisation
     # Lower, central, and upper values for each parameter
@@ -132,12 +138,12 @@ def run_opt_sim(args, numbermodelruns, inputfile, usernamespace):
             break
 
         # Stop optimisation if successive fitness values are within a percentage threshold
-#        if iteration > 2:
-#            fitnessvaluesclose = (np.abs(fitnessvalueshist[iteration - 2] - fitnessvalueshist[iteration - 1]) / fitnessvalueshist[iteration - 1]) * 100
-#            fitnessvaluesthres = 0.1
-#            if fitnessvaluesclose < fitnessvaluesthres:
-#                print('\nTaguchi optimisation stopped as successive fitness values within {}%'.format(fitnessvaluesthres))
-#                break
+        if iteration > 2:
+            fitnessvaluesclose = (np.abs(fitnessvalueshist[iteration - 2] - fitnessvalueshist[iteration - 1]) / fitnessvalueshist[iteration - 1]) * 100
+            fitnessvaluesthres = 0.1
+            if fitnessvaluesclose < fitnessvaluesthres:
+                print('\nTaguchi optimisation stopped as successive fitness values within {}%'.format(fitnessvaluesthres))
+                break
 
     # Save optimisation parameters history and fitness values history to file
     opthistfile = inputfileparts[0] + '_hist.pickle'
@@ -164,6 +170,9 @@ def taguchi_code_blocks(inputfile, taguchinamespace):
         # Strip out any newline characters and comments that must begin with double hashes
         inputlines = [line.rstrip() for line in f if(not line.startswith('##') and line.rstrip('\n'))]
     
+    # Store length of dict
+    taglength = len(taguchinamespace)
+    
     x = 0
     while(x < len(inputlines)):
         if(inputlines[x].startswith('#taguchi:')):
@@ -184,6 +193,10 @@ def taguchi_code_blocks(inputfile, taguchinamespace):
             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
 
@@ -202,6 +215,9 @@ def construct_OA(optparams):
         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
@@ -215,7 +231,7 @@ def construct_OA(optparams):
     
     # Load the appropriate OA
     if k <= 4:
-        OA = np.load(os.path.join('user_libs', 'OA_9_4_3_2.npy'))
+        OA = np.load(os.path.join(oadirectory, 'OA_9_4_3_2.npy'))
 
         # Number of experiments
         N = OA.shape[0]
@@ -227,7 +243,7 @@ def construct_OA(optparams):
         OA = OA[:, 0:k]
 
     elif k <= 7:
-        OA = np.load(os.path.join('user_libs',  'OA_18_7_3_2.npy'))
+        OA = np.load(os.path.join(oadirectory, 'OA_18_7_3_2.npy'))
 
         # Number of experiments
         N = OA.shape[0]
@@ -238,8 +254,10 @@ def construct_OA(optparams):
         # Cut down OA columns to number of parameters to optimise
         OA = OA[:, 0:k]
 
-    # THIS CASE NEEDS FURTHER TESTING
     else:
+        # THIS CASE NEEDS FURTHER TESTING
+        print('\nTaguchi optimisation, WARNING: Optimising more than 7 parameters is currently an experimental feature!')
+              
         p = int(np.ceil(np.log(k * (s - 1) + 1) / np.log(s)))
         
         # Number of experiments
@@ -380,8 +398,8 @@ def calculate_optimal_levels(optparams, levels, levelsopt, fitnessvalues, OA, N,
         # Calculate optimal level from table of responses
         optlevel = np.where(responses == np.amax(responses))[0]
         
-        # If 2 experiments produce the same fitness value (this shouldn't happen if the fitness function is designed correctly)
-        if len(optlevel):
+        # If 2 experiments produce the same fitness value (this shouldn't happen if the fitness function is designed correctly) pick first level
+        if len(optlevel) > 1:
             optlevel = optlevel[0]
         
         levelsopt[p] = optlevel

gprMax/pml.py

@@ -23,7 +23,7 @@ from gprMax.pml_1order_update import *
 from gprMax.pml_2order_update import *
 
 
-class CFSParameter:
+class CFSParameter(object):
     """Individual CFS parameter (e.g. alpha, kappa, or sigma)."""
     
     # Allowable scaling profiles and directions
@@ -49,7 +49,7 @@ class CFSParameter:
         self.max = max
 
 
-class CFS:
+class CFS(object):
     """CFS term for PML."""
     
     def __init__(self):
@@ -140,7 +140,7 @@ class CFS:
         return Evalues, Hvalues
 
 
-class PML:
+class PML(object):
     """PML - the implementation comes from the derivation in: http://dx.doi.org/10.1109/TAP.2011.2180344"""
     
     directions = {0: 'xminus', 1: 'yminus', 2: 'zminus', 3: 'xplus', 4: 'yplus', 5: 'zplus'}
@@ -165,6 +165,10 @@ class PML:
         self.CFS = G.cfs
         if not self.CFS:
             self.CFS = [CFS()]
+        self.initialise_field_arrays()
+        
+    def initialise_field_arrays(self):
+        """Initialise arrays to store fields in PML."""
         
         # Subscript notation, e.g. 'EPhiyxz' means the electric field Phi vector, of which the
         # component being corrected is y, the stretching direction is x, and field derivative
@@ -188,15 +192,6 @@ class PML:
             self.HPhixzy = np.zeros((len(self.CFS), self.nx + 1, self.ny, self.nz), dtype=floattype)
             self.HPhiyzx = np.zeros((len(self.CFS), self.nx, self.ny + 1, self.nz), dtype=floattype)
 
-        self.ERA = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
-        self.ERB = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
-        self.ERE = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
-        self.ERF = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
-        self.HRA = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
-        self.HRB = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
-        self.HRE = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
-        self.HRF = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
-                
     def calculate_update_coeffs(self, er, mr, G):
         """Calculates electric and magnetic update coefficients for the PML.
             
@@ -206,6 +201,15 @@ class PML:
             G (class): Grid class instance - holds essential parameters describing the model.
         """
         
+        self.ERA = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
+        self.ERB = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
+        self.ERE = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
+        self.ERF = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
+        self.HRA = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
+        self.HRB = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
+        self.HRE = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
+        self.HRF = np.zeros((len(self.CFS), self.thickness + 1), dtype=floattype)
+        
         for x, cfs in enumerate(self.CFS):
             if not cfs.sigma.max:
                 cfs.calculate_sigmamax(self.direction, er, mr, G)

gprMax/receivers.py

@@ -16,21 +16,18 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
-class Rx:
+class Rx(object):
     """Receiever output points."""
     
     availableoutputs = ['Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', 'Ix', 'Iy', 'Iz']
     
-    def __init__(self, xcoord=None, ycoord=None, zcoord=None):
-        """
-        Args:
-            xcoord (float): x-coordinate of location in model.
-            ycoord (float): y-coordinate of location in model.
-            zcoord (float): z-coordinate of location in model.
-        """
+    def __init__(self):
         
         self.ID = None
         self.outputs = []
-        self.xcoord = xcoord
-        self.ycoord = ycoord
-        self.zcoord = zcoord
+        self.xcoord = None
+        self.ycoord = None
+        self.zcoord = None
+        self.xcoordbase = None
+        self.ycoordbase = None
+        self.zcoordbase = None

gprMax/snapshots.py

@@ -16,7 +16,7 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
-import sys
+import os, sys
 import numpy as np
 from struct import pack
 
@@ -25,7 +25,7 @@ from gprMax.grid import Ix, Iy, Iz
 from gprMax.utilities import round_value
 
 
-class Snapshot:
+class Snapshot(object):
     """Snapshots of the electric and magnetic field values."""
     
     # Set string for byte order
@@ -61,7 +61,7 @@ class Snapshot:
         self.dy = dy
         self.dz = dz
         self.time = time
-        self.filename = filename
+        self.basefilename = filename
 
     def prepare_vtk_imagedata(self, modelrun, numbermodelruns, G):
         """Prepares a VTK ImageData (.vti) file for a snapshot.
@@ -77,11 +77,15 @@ class Snapshot:
         self.vtk_ny = self.yf - self.ys
         self.vtk_nz = self.zf - self.zs
         
-        # Construct filename from user-supplied name and model run number
+        # Create directory and construct filename from user-supplied name and model run number
         if numbermodelruns == 1:
-            self.filename = G.inputdirectory + self.filename + '.vti'
+            snapshotdir = os.path.join(G.inputdirectory, os.path.splitext(G.inputfilename)[0] + '_snaps')
         else:
-            self.filename = G.inputdirectory + self.filename + '_' + str(modelrun) + '.vti'
+            snapshotdir = os.path.join(G.inputdirectory, os.path.splitext(G.inputfilename)[0] + '_snaps' + str(modelrun))
+        
+        if not os.path.exists(snapshotdir):
+            os.mkdir(snapshotdir)
+        self.filename = os.path.abspath(os.path.join(snapshotdir, self.basefilename + '.vti'))
         
         # Calculate number of cells according to requested sampling
         self.vtk_xscells = round_value(self.xs / self.dx)

gprMax/sources.py

@@ -24,20 +24,30 @@ from gprMax.grid import Ix, Iy, Iz
 from gprMax.utilities import round_value
 
 
-class VoltageSource:
-    """The voltage source can be a hard source if it's resistance is zero, i.e. the time variation of the specified electric field component is prescribed. If it's resistance is non-zero it behaves as a resistive voltage source."""
-    
+class Source(object):
+    """Super-class which describes a generic source."""
+
     def __init__(self):
         self.ID = None
         self.polarisation = None
         self.xcoord = None
         self.ycoord = None
         self.zcoord = None
+        self.xcoordbase = None
+        self.ycoordbase = None
+        self.zcoordbase = None
         self.start = None
         self.stop = None
-        self.resistance = None
         self.waveformID = None
 
+
+class VoltageSource(Source):
+    """The voltage source can be a hard source if it's resistance is zero, i.e. the time variation of the specified electric field component is prescribed. If it's resistance is non-zero it behaves as a resistive voltage source."""
+    
+    def __init__(self):
+        super(Source, self).__init__()
+        self.resistance = None
+
     def update_electric(self, abstime, updatecoeffsE, ID, Ex, Ey, Ez, G):
         """Updates electric field values for a voltage source.
             
@@ -59,19 +69,22 @@ class VoltageSource:
             
             if self.polarisation is 'x':
                 if self.resistance != 0:
-                    Ex[i, j, k] -= updatecoeffsE[ID[0, i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (self.resistance * G.dy * G.dz))
+                    componentID = 'E' + self.polarisation
+                    Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (self.resistance * G.dy * G.dz))
                 else:
                     Ex[i, j, k] = -1 * waveform.amp * waveform.calculate_value(time, G.dt) / G.dx
 
             elif self.polarisation is 'y':
                 if self.resistance != 0:
-                    Ey[i, j, k] -= updatecoeffsE[ID[1, i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (self.resistance * G.dx * G.dz))
+                    componentID = 'E' + self.polarisation
+                    Ey[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (self.resistance * G.dx * G.dz))
                 else:
                     Ey[i, j, k] = -1 * waveform.amp * waveform.calculate_value(time, G.dt) / G.dy
 
             elif self.polarisation is 'z':
                 if self.resistance != 0:
-                    Ez[i, j, k] -= updatecoeffsE[ID[2, i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (self.resistance * G.dx * G.dy))
+                    componentID = 'E' + self.polarisation
+                    Ez[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (self.resistance * G.dx * G.dy))
                 else:
                     Ez[i, j, k] = -1 * waveform.amp * waveform.calculate_value(time, G.dt) / G.dz
 
@@ -83,8 +96,12 @@ class VoltageSource:
         """
 
         if self.resistance != 0:
-            ID = 'E' + self.polarisation
-            requirednumID = G.ID[G.IDlookup[ID], self.xcoord, self.ycoord, self.zcoord]
+            i = self.xcoord
+            j = self.ycoord
+            k = self.zcoord
+            
+            componentID = 'E' + self.polarisation
+            requirednumID = G.ID[G.IDlookup[componentID], i, j, k]
             material = next(x for x in G.materials if x.numID == requirednumID)
             newmaterial = deepcopy(material)
             newmaterial.ID = material.ID + '+VoltageSource_' + str(self.resistance)
@@ -99,22 +116,15 @@ class VoltageSource:
             elif self.polarisation == 'z':
                 newmaterial.se += G.dz / (self.resistance * G.dx * G.dy)
 
-            G.ID[G.IDlookup[ID], self.xcoord, self.ycoord, self.zcoord] = newmaterial.numID
+            G.ID[G.IDlookup[componentID], i, j, k] = newmaterial.numID
             G.materials.append(newmaterial)
 
 
-class HertzianDipole:
+class HertzianDipole(Source):
     """The Hertzian dipole is an additive source (electric current density)."""
     
     def __init__(self):
-        self.ID = None
-        self.polarisation = None
-        self.xcoord = None
-        self.ycoord = None
-        self.zcoord = None
-        self.start = None
-        self.stop = None
-        self.waveformID = None
+        super(Source, self).__init__()
 
     def update_electric(self, abstime, updatecoeffsE, ID, Ex, Ey, Ez, G):
         """Updates electric field values for a Hertzian dipole.
@@ -136,27 +146,23 @@ class HertzianDipole:
             waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
             
             if self.polarisation is 'x':
-                Ex[i, j, k] -= updatecoeffsE[ID[0, i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dy * G.dz))
+                componentID = 'E' + self.polarisation
+                Ex[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dy * G.dz))
 
             elif self.polarisation is 'y':
-                Ey[i, j, k] -= updatecoeffsE[ID[1, i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dx * G.dz))
+                componentID = 'E' + self.polarisation
+                Ey[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dx * G.dz))
 
             elif self.polarisation is 'z':
-                Ez[i, j, k] -= updatecoeffsE[ID[2, i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dx * G.dy))
+                componentID = 'E' + self.polarisation
+                Ez[i, j, k] -= updatecoeffsE[ID[G.IDlookup[componentID], i, j, k], 4] * waveform.amp * waveform.calculate_value(time, G.dt) * (1 / (G.dx * G.dy))
 
 
-class MagneticDipole:
+class MagneticDipole(Source):
     """The magnetic dipole is an additive source (magnetic current density)."""
     
     def __init__(self):
-        self.ID = None
-        self.polarisation = None
-        self.xcoord = None
-        self.ycoord = None
-        self.zcoord = None
-        self.start = None
-        self.stop = None
-        self.waveformID = None
+        super(Source, self).__init__()
 
     def update_magnetic(self, abstime, updatecoeffsH, ID, Hx, Hy, Hz, G):
         """Updates magnetic field values for a magnetic dipole.
@@ -187,7 +193,7 @@ class MagneticDipole:
                 Hz[i, j, k] -= waveform.amp  * waveform.calculate_value(time, G.dt) * (G.dt / (G.dx * G.dy * G.dz))
 
 
-class TransmissionLine:
+class TransmissionLine(Source):
     """The transmission line source is a one-dimensional transmission line which is attached virtually to a grid cell."""
     
     def __init__(self, G):
@@ -196,15 +202,8 @@ class TransmissionLine:
             G (class): Grid class instance - holds essential parameters describing the model.
         """
         
-        self.ID = None
-        self.polarisation = None
-        self.xcoord = None
-        self.ycoord = None
-        self.zcoord = None
-        self.start = None
-        self.stop = None
+        super(Source, self).__init__()
         self.resistance = None
-        self.waveformID = None
         
         # Coefficients for ABC termination of end of the transmission line
         self.abcv0 = 0

gprMax/utilities.py

@@ -21,7 +21,7 @@ import decimal as d
 from pyfiglet import Figlet
 
 
-class ListStream:
+class ListStream(object):
     """A list can be streamed into. Required when temporarily redirecting stdio to capture output from users Python code blocks."""
     
     def __init__(self):

gprMax/waveforms.py

@@ -21,7 +21,7 @@ import numpy as np
 from gprMax.utilities import round_value
 
 
-class Waveform:
+class Waveform(object):
     """Definitions of waveform shapes that can be used with sources."""
     
     types = ['gaussian', 'gaussiandot', 'gaussiandotnorm', 'gaussiandotdot', 'gaussiandotdotnorm', 'ricker', 'sine', 'contsine', 'impulse', 'user']
@@ -89,7 +89,7 @@ class Waveform:
             waveform = ramp * np.sin(2 * np.pi * self.freq * time)
 
         elif self.type == 'impulse':
-            # time < G.dt condition required to do impulsive magnetic dipole
+            # time < dt condition required to do impulsive magnetic dipole
             if time == 0 or time < dt:
                 waveform = 1
             elif time >= dt:

gprMax/writer_hdf5.py

@@ -19,6 +19,7 @@
 import h5py
 import numpy as np
 
+import gprMax
 from gprMax.constants import floattype
 from gprMax.grid import Ix, Iy, Iz
 
@@ -35,6 +36,7 @@ def prepare_hdf5(outputfile, G):
     """
 
     f = h5py.File(outputfile, 'w')
+    f.attrs['gprMax'] = gprMax.__version__
     f.attrs['Title'] = G.title
     f.attrs['Iterations'] = G.iterations
     f.attrs['nx, ny, nz'] = (G.nx, G.ny, G.nz)

gprMax/yee_cell_build.pyx

@@ -108,10 +108,10 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
     """
     
     cdef Py_ssize_t i, j, k
-    cdef int numID1, numID2, numID3, numID4
+    cdef int numID1, numID2, numID3, numID4, componentID
 
     # Ex component
-    componentID = 0
+    componentID = G.IDlookup['Ex']
     for i in range(0, G.nx):
         for j in range(1, G.ny):
             for k in range(1, G.nz):
@@ -133,7 +133,7 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                         create_electric_average(i, j, k, numID1, numID2, numID3, numID4, componentID, G)
 
     # Ey component
-    componentID = 1
+    componentID = G.IDlookup['Ey']
     for i in range(1, G.nx):
         for j in range(0, G.ny):
             for k in range(1, G.nz):
@@ -155,7 +155,7 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                         create_electric_average(i, j, k, numID1, numID2, numID3, numID4, componentID, G)
 
     # Ez component
-    componentID = 2
+    componentID = G.IDlookup['Ez']
     for i in range(1, G.nx):
         for j in range(1, G.ny):
             for k in range(0, G.nz):
@@ -185,10 +185,10 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
     """
     
     cdef Py_ssize_t i, j, k
-    cdef int numID1, numID2
+    cdef int numID1, numID2, componentID
 
     # Hx component
-    componentID = 3
+    componentID = G.IDlookup['Hx']
     for i in range(1, G.nx):
         for j in range(0, G.ny):
             for k in range(0, G.nz):
@@ -208,7 +208,7 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                         create_magnetic_average(i, j, k, numID1, numID2, componentID, G)
 
     # Hy component
-    componentID = 4
+    componentID = G.IDlookup['Hy']
     for i in range(0, G.nx):
         for j in range(1, G.ny):
             for k in range(0, G.nz):
@@ -228,7 +228,7 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                         create_magnetic_average(i, j, k, numID1, numID2, componentID, G)
 
     # Hz component
-    componentID = 5
+    componentID = G.IDlookup['Hz']
     for i in range(0, G.nx):
         for j in range(0, G.ny):
             for k in range(1, G.nz):