Work on output directory, snapshots, and geometry view path settings.

gprMax/cmds_multiple.py

@@ -664,7 +664,7 @@ class Snapshot(UserObjectMulti):
         #else:
         s = SnapshotUser(xs, ys, zs, xf, yf, zf, dx, dy, dz, iterations, filename)
 
-        log.info(f'Snapshot from {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m, discretisation {dx * grid.dx:g}m, {dy * grid.dy:g}m, {dz * grid.dz:g}m, at {s.time * grid.dt:g} secs with filename {s.basefilename} created.')
+        log.info(f'Snapshot from {xs * grid.dx:g}m, {ys * grid.dy:g}m, {zs * grid.dz:g}m, to {xf * grid.dx:g}m, {yf * grid.dy:g}m, {zf * grid.dz:g}m, discretisation {dx * grid.dx:g}m, {dy * grid.dy:g}m, {dz * grid.dz:g}m, at {s.time * grid.dt:g} secs with filename {s.filename} created.')
 
         grid.snapshots.append(s)
 

gprMax/cmds_single_use.py

@@ -504,7 +504,7 @@ class OutputDir(UserObjectSingle):
         self.order = 11
 
     def create(self, grid, uip):
-        grid.outputdirectory = self.kwargs['dir']
+        config.model_configs[grid.model_num].set_output_file_path(self.kwargs['dir'])
 
 
 class NumberOfModelRuns(UserObjectSingle):

gprMax/config.py

@@ -73,22 +73,20 @@ class ModelConfig:
 
         self.reuse_geometry = False
 
+        # String to print at start of each model run
+        inputfilestr = f'\n--- Model {self.i + 1}/{sim_config.model_end}, input file: {sim_config.input_file_path}'
+        self.set_inputfilestr(inputfilestr)
+
         if not sim_config.single_model:
             self.appendmodelnumber = str(self.i + 1) # Indexed from 1
         else:
             self.appendmodelnumber = ''
 
         # Output file path for specific model
-        parts = sim_config.output_file_path.with_suffix('').parts
+        self.set_output_file_path()
-        self.output_file_path = Path(*parts[:-1], parts[-1] + self.appendmodelnumber)
-        self.output_file_path_ext = self.output_file_path.with_suffix('.out')
 
-        # Make a snapshot directory
+        # Specify a snapshot directory
-        self.snapshot_dir = '_snaps'
+        self.set_snapshots_file_path()
-
-        # String to print at start of each model run
-        inputfilestr = f'\n--- Model {self.i + 1}/{sim_config.model_end}, input file: {sim_config.input_file_path}'
-        self.set_inputfilestr(inputfilestr)
 
         # Numerical dispersion analysis parameters
         #   highestfreqthres: threshold (dB) down from maximum power (0dB) of main frequency used
@@ -129,6 +127,34 @@ class ModelConfig:
         """
         self.inputfilestr = Fore.GREEN + f"{inputfilestr} {'-' * (get_terminal_width() - 1 - len(inputfilestr))}\n" + Style.RESET_ALL
 
+    def set_output_file_path(self, outputdir=None):
+        """Output file path can be provided by the user via the API or an input file
+            command. If they haven't provided one use the input file path instead.
+
+        Args:
+            outputdir (str): Output file directory given from input file command.
+        """
+
+        if not outputdir:
+            try:
+                self.output_file_path = Path(self.args.outputfile)
+            except AttributeError:
+                self.output_file_path = sim_config.input_file_path.with_suffix('')
+        else:
+            try:
+                Path(outputdir).mkdir(exist_ok=True)
+                self.output_file_path = Path(outputdir, sim_config.input_file_path.stem)
+            except AttributeError:
+                self.output_file_path = sim_config.input_file_path.with_suffix('')
+
+        parts = self.output_file_path.parts
+        self.output_file_path = Path(*parts[:-1], parts[-1] + self.appendmodelnumber)
+        self.output_file_path_ext = self.output_file_path.with_suffix('.out')
+
+    def set_snapshots_file_path(self):
+        """Set directory to store any snapshots."""
+        parts = self.output_file_path.with_suffix('').parts
+        self.snapshot_file_path = Path(*parts[:-1], parts[-1] + '_snaps')
 
 class SimulationConfig:
     """Configuration parameters for a standard simulation.
@@ -200,8 +226,8 @@ class SimulationConfig:
 
         # Set more complex parameters
         self.set_precision()
+        self.get_byteorder()
         self.set_input_file_path()
-        self.set_output_file_path()
         self.set_model_start_end()
         self.set_single_model()
 
@@ -237,14 +263,22 @@ class SimulationConfig:
                       'cython_float_or_double': cython.float,
                       'cython_complex': cython.floatcomplex,
                       'C_float_or_double': 'float',
-                      'C_complex': 'pycuda::complex<float>'}
+                      'C_complex': 'pycuda::complex<float>',
+                      'vtk_float': 'Float32'}
         elif self.general['precision'] == 'double':
             self.dtypes = {'float_or_double': np.float64,
                       'complex': np.complex128,
                       'cython_float_or_double': cython.double,
                       'cython_complex': cython.doublecomplex,
                       'C_float_or_double': 'double',
-                      'C_complex': 'pycuda::complex<double>'}
+                      'C_complex': 'pycuda::complex<double>',
+                      'vtk_float': 'Float64'}
+
+    def get_byteorder(self):
+        """Check the byte order of system to use for VTK files, i.e. geometry
+            views and snapshots.
+        """
+        self.vtk_byteorder = 'LittleEndian' if sys.byteorder == 'little' else 'BigEndian'
 
     def set_single_model(self):
         if self.model_start == 0 and self.model_end == 1:
@@ -269,20 +303,12 @@ class SimulationConfig:
         self.model_end = modelend
 
     def set_input_file_path(self):
-        """If the API is in use an id for the simulation must be provided."""
+        """Set input file path for CLI or API."""
         if self.args.inputfile is None:
             self.input_file_path = Path(self.args.outputfile)
         else:
             self.input_file_path = Path(self.args.inputfile)
 
-    def set_output_file_path(self):
-        """Output file path can be provided by the user. If they havent provided one
-            use the inputfile file path instead."""
-        try:
-            self.output_file_path = Path(self.args.outputfile)
-        except AttributeError:
-            self.output_file_path = Path(self.input_file_path)
-
 
 class SimulationConfigMPI(SimulationConfig):
     """Configuration parameters for a MPI simulation.

gprMax/fields_outputs.py

@@ -31,7 +31,7 @@ def store_outputs(G):
         iteration (int): Current iteration number.
         Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic
                                                 field values.
-        G (FDTDGrid): Holds essential parameters describing a model.
+        G (FDTDGrid): Parameters describing a grid in a model.
     """
 
     iteration = G.iteration
@@ -105,7 +105,7 @@ def write_hdf5_main_grid_outputfile(outputfile, G):
 
     Args:
         outputfile (str): Name of the output file.
-        G (FDTDGrid): Holds essential parameters describing a model.
+        G (FDTDGrid): Parameters describing a grid in a model.
     """
 
     write_data(outputfile, G)
@@ -116,7 +116,7 @@ def write_hdf5_sub_grid_outputfile(outputfile, G):
 
     Args:
         outputfile (str): Name of the output file.
-        G (FDTDGrid): Holds essential parameters describing a model.
+        G (FDTDGrid): Parameters describing a grid in a model.
     """
 
     stem = outputfile.stem
@@ -145,7 +145,7 @@ def write_data(outputfile, G):
 
     Args:
         outputfile (str): Name of the output file.
-        G (FDTDGrid): Holds essential parameters describing a model.
+        G (FDTDGrid): Parameters describing a grid in a model.
     """
 
     f = h5py.File(outputfile, 'w')

gprMax/geometry_outputs.py

@@ -16,6 +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 logging
 import os
 from pathlib import Path
 import sys
@@ -31,16 +32,13 @@ from .cython.geometry_outputs import define_normal_geometry
 from .cython.geometry_outputs import define_fine_geometry
 from .utilities import round_value
 
+log = logging.getLogger(__name__)
+
 
 class GeometryView:
     """Views of the geometry of the model."""
 
-    if sys.byteorder == 'little':
+    def __init__(self, xs=None, ys=None, zs=None, xf=None, yf=None, zf=None, dx=None, dy=None, dz=None, filename=None, fileext=None, G=None):
-        byteorder = 'LittleEndian'
-    else:
-        byteorder = 'BigEndian'
-
-    def __init__(self, xs=None, ys=None, zs=None, xf=None, yf=None, zf=None, dx=None, dy=None, dz=None, filename=None, fileext=None, grid=None):
         """
         Args:
             xs, xf, ys, yf, zs, zf (int): Extent of the volume in cells.
@@ -48,6 +46,7 @@ class GeometryView:
             filename (str): Filename to save to.
             fileext (str): File extension of VTK file - either '.vti' for a per cell
                     geometry view, or '.vtp' for a per cell edge geometry view.
+            G (FDTDGrid): Parameters describing a grid in a model.
         """
 
         self.xs = xs
@@ -64,6 +63,7 @@ class GeometryView:
         self.dz = dz
         self.basefilename = filename
         self.fileext = fileext
+        self.G = G
 
         if self.fileext == '.vti':
             # Calculate number of cells according to requested sampling for geometry view
@@ -77,9 +77,9 @@ class GeometryView:
             self.vtk_nycells = round_value(self.ny / self.dy)
             self.vtk_nzcells = round_value(self.nz / self.dz)
             self.vtk_ncells = self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells
-            self.datawritesize = (np.dtype(np.uint32).itemsize * self.vtk_ncells
+            self.datawritesize = (np.dtype(np.uint32).itemsize * self.vtk_ncells +
-                                  + 2 * np.dtype(np.int8).itemsize * self.vtk_ncells
+                                  2 * np.dtype(np.int8).itemsize * self.vtk_ncells +
-                                  + 3 * np.dtype(np.uint32).itemsize)
+                                  3 * np.dtype(np.uint32).itemsize)
 
         elif self.fileext == '.vtp':
             self.vtk_numpoints = (self.nx + 1) * (self.ny + 1) * (self.nz + 1)
@@ -89,45 +89,44 @@ class GeometryView:
             self.vtk_nylines = self.ny * (self.nx + 1) * (self.nz + 1)
             self.vtk_nzlines = self.nz * (self.nx + 1) * (self.ny + 1)
             self.vtk_numlines = self.vtk_nxlines + self.vtk_nylines + self.vtk_nzlines
-            self.vtk_connectivity_offset = round_value((self.vtk_numpoints
+            self.vtk_connectivity_offset = round_value((self.vtk_numpoints *
-                                                        * self.vtk_numpoint_components
+                                                       self.vtk_numpoint_components *
-                                                        * np.dtype(np.float32).itemsize)
+                                                       np.dtype(np.float32).itemsize) +
-                                                       + np.dtype(np.uint32).itemsize)
+                                                       np.dtype(np.uint32).itemsize)
-            self.vtk_offsets_offset = round_value(self.vtk_connectivity_offset
+            self.vtk_offsets_offset = round_value(self.vtk_connectivity_offset +
-                                                  + (self.vtk_numlines * self.vtk_numline_components * np.dtype(np.uint32).itemsize)
+                                                 (self.vtk_numlines *
-                                                  + np.dtype(np.uint32).itemsize)
+                                                 self.vtk_numline_components *
-            self.vtk_materials_offset = round_value(self.vtk_offsets_offset
+                                                 np.dtype(np.uint32).itemsize) +
-                                                    + (self.vtk_numlines * np.dtype(np.uint32).itemsize)
+                                                 np.dtype(np.uint32).itemsize)
-                                                    + np.dtype(np.uint32).itemsize)
+            self.vtk_materials_offset = round_value(self.vtk_offsets_offset +
-            vtk_cell_offsets = ((self.vtk_numline_components * self.vtk_numlines)
+                                                   (self.vtk_numlines *
-                                + self.vtk_numline_components - self.vtk_numline_components - 1) // self.vtk_numline_components + 1
+                                                   np.dtype(np.uint32).itemsize) +
-            self.datawritesize = (np.dtype(np.float32).itemsize * self.vtk_numpoints * self.vtk_numpoint_components
+                                                   np.dtype(np.uint32).itemsize)
-                                  + np.dtype(np.uint32).itemsize * self.vtk_numlines * self.vtk_numline_components
+            vtk_cell_offsets = (((self.vtk_numline_components * self.vtk_numlines) +
-                                  + np.dtype(np.uint32).itemsize * self.vtk_numlines
+                               self.vtk_numline_components - self.vtk_numline_components - 1) //
-                                  + np.dtype(np.uint32).itemsize * vtk_cell_offsets
+                               self.vtk_numline_components + 1)
-                                  + np.dtype(np.uint32).itemsize * 4)
+            self.datawritesize = (np.dtype(np.float32).itemsize * self.vtk_numpoints *
+                                 self.vtk_numpoint_components + np.dtype(np.uint32).itemsize *
+                                 self.vtk_numlines * self.vtk_numline_components +
+                                 np.dtype(np.uint32).itemsize * self.vtk_numlines +
+                                 np.dtype(np.uint32).itemsize * vtk_cell_offsets +
+                                 np.dtype(np.uint32).itemsize * 4)
 
-    def set_filename(self, appendmodelnumber):
+    def set_filename(self):
-        """
+        """Construct filename from user-supplied name and model run number."""
-        Construct filename from user-supplied name and model run number.
+        parts = config.model_configs[self.G.model_num].output_file_path.with_suffix('').parts
-
+        self.filename = Path(*parts[:-1], parts[-1])
-        Args:
-            appendmodelnumber (str): Text to append to filename.
-        """
-
-        parts = config.sim_config.input_file_path.parts
-        self.filename = Path(*parts[:-1], parts[-1] + appendmodelnumber)
         self.filename = self.filename.with_suffix(self.fileext)
 
     def write_vtk(self, G, pbar):
-        """
+        """Writes the geometry information to a VTK file.
-        Writes the geometry information to a VTK file. Either ImageData (.vti) for a
+            Either ImageData (.vti) for a per-cell geometry view, or
-        per-cell geometry view, or PolygonalData (.vtp) for a per-cell-edge geometry view.
+            PolygonalData (.vtp) for a per-cell-edge geometry view.
 
             N.B. No Python 3 support for VTK at time of writing (03/2015)
 
         Args:
-            G (class): Grid class instance - holds essential parameters describing the model.
+            G (FDTDGrid): Parameters describing a grid in a model.
             pbar (class): Progress bar class instance.
         """
 
@@ -143,18 +142,29 @@ class GeometryView:
             for index, rx in enumerate(G.rxs):
                 self.rxs[rx.xcoord, rx.ycoord, rx.zcoord] = index + 1
 
-            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_srcs_pml_offset = round_value((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)
+                                               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'))
-                f.write('<VTKFile type="ImageData" version="1.0" byte_order="{}">\n'.format(GeometryView.byteorder).encode('utf-8'))
+                f.write(f'<VTKFile type="ImageData" version="1.0" byte_order="{config.sim_config.vtk_byteorder}">\n'.encode('utf-8'))
-                f.write('<ImageData WholeExtent="{} {} {} {} {} {}" Origin="0 0 0" Spacing="{:.3} {:.3} {:.3}">\n'.format(self.vtk_xscells, self.vtk_xfcells, self.vtk_yscells, self.vtk_yfcells, self.vtk_zscells, self.vtk_zfcells, self.dx * G.dx, self.dy * G.dy, self.dz * G.dz).encode('utf-8'))
+                f.write(f'<ImageData WholeExtent="{self.vtk_xscells} {self.vtk_xfcells} {self.vtk_yscells} {self.vtk_yfcells} {self.vtk_zscells} {self.vtk_zfcells}" Origin="0 0 0" Spacing="{self.dx * G.dx:.3} {self.dy * G.dy:.3} {self.dz * G.dz:.3}">\n'.encode('utf-8'))
-                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(f'<Piece Extent="{self.vtk_xscells} {self.vtk_xfcells} {self.vtk_yscells} {self.vtk_yfcells} {self.vtk_zscells} {self.vtk_zfcells}">\n'.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_PML" format="appended" offset="{}" />\n'.format(vtk_srcs_pml_offset).encode('utf-8'))
+                f.write(f'<DataArray type="Int8" Name="Sources_PML" format="appended" offset="{vtk_srcs_pml_offset}" />\n'.encode('utf-8'))
-                f.write('<DataArray type="Int8" Name="Receivers" format="appended" offset="{}" />\n'.format(vtk_rxs_offset).encode('utf-8'))
+                f.write(f'<DataArray type="Int8" Name="Receivers" format="appended" offset="{vtk_rxs_offset}" />\n'.encode('utf-8'))
                 f.write('</CellData>\n'.encode('utf-8'))
                 f.write('</Piece>\n</ImageData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
 
@@ -207,33 +217,30 @@ class GeometryView:
         elif self.fileext == '.vtp':
             with open(self.filename, 'wb') as f:
                 f.write('<?xml version="1.0"?>\n'.encode('utf-8'))
-                f.write('<VTKFile type="PolyData" version="1.0" byte_order="{}">\n'.format(GeometryView.byteorder).encode('utf-8'))
+                f.write(f'<VTKFile type="PolyData" version="1.0" byte_order="{config.sim_config.vtk_byteorder}">\n'.encode('utf-8'))
-                f.write('<PolyData>\n<Piece NumberOfPoints="{}" NumberOfVerts="0" NumberOfLines="{}" NumberOfStrips="0" NumberOfPolys="0">\n'.format(self.vtk_numpoints, self.vtk_numlines).encode('utf-8'))
+                f.write(f'<PolyData>\n<Piece NumberOfPoints="{self.vtk_numpoints}" NumberOfVerts="0" NumberOfLines="{self.vtk_numlines}" NumberOfStrips="0" NumberOfPolys="0">\n'.encode('utf-8'))
-
                 f.write('<Points>\n<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="0" />\n</Points>\n'.encode('utf-8'))
-                f.write('<Lines>\n<DataArray type="UInt32" Name="connectivity" format="appended" offset="{}" />\n'.format(self.vtk_connectivity_offset).encode('utf-8'))
+                f.write(f'<Lines>\n<DataArray type="UInt32" Name="connectivity" format="appended" offset="{self.vtk_connectivity_offset}" />\n'.encode('utf-8'))
-                f.write('<DataArray type="UInt32" Name="offsets" format="appended" offset="{}" />\n</Lines>\n'.format(self.vtk_offsets_offset).encode('utf-8'))
+                f.write(f'<DataArray type="UInt32" Name="offsets" format="appended" offset="{self.vtk_offsets_offset}" />\n</Lines>\n'.encode('utf-8'))
-
                 f.write('<CellData Scalars="Material">\n'.encode('utf-8'))
                 f.write('<DataArray type="UInt32" Name="Material" format="appended" offset="{}" />\n'.format(self.vtk_materials_offset).encode('utf-8'))
                 f.write('</CellData>\n'.encode('utf-8'))
-
                 f.write('</Piece>\n</PolyData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
 
                 # Coordinates of each point
                 points = np.zeros((self.vtk_numpoints, 3), dtype=np.float32)
 
-                # Number of x components
-
                 # Node connectivity. Each index contains a pair of connected x nodes
                 x_lines = np.zeros((self.vtk_nxlines, 2), dtype=np.uint32)
                 # Material at Ex location in Yee cell.
                 x_materials = np.zeros((self.vtk_nxlines), dtype=np.uint32)
-
+                # Node connectivity. Each index contains a pair of connected y nodes
                 y_lines = np.zeros((self.vtk_nylines, 2), dtype=np.uint32)
+                # Material at Ey location in Yee cell.
                 y_materials = np.zeros((self.vtk_nylines), dtype=np.uint32)
-
+                # Node connectivity. Each index contains a pair of connected z nodes
                 z_lines = np.zeros((self.vtk_nzlines, 2), dtype=np.uint32)
+                # Material at Ez location in Yee cell.
                 z_materials = np.zeros((self.vtk_nzlines), dtype=np.uint32)
 
                 define_fine_geometry(self.nx,
@@ -263,7 +270,8 @@ class GeometryView:
                 pbar.update(n=points.nbytes)
 
                 # Write connectivity data
-                f.write(pack('I', np.dtype(np.uint32).itemsize * self.vtk_numlines * self.vtk_numline_components))
+                f.write(pack('I', np.dtype(np.uint32).itemsize *
+                                  self.vtk_numlines * self.vtk_numline_components))
                 pbar.update(n=4)
                 f.write(x_lines)
                 pbar.update(n=x_lines.nbytes)
@@ -275,7 +283,9 @@ class GeometryView:
                 # Write cell type (line) offsets
                 f.write(pack('I', np.dtype(np.uint32).itemsize * self.vtk_numlines))
                 pbar.update(n=4)
-                for vtk_offsets in range(self.vtk_numline_components, (self.vtk_numline_components * self.vtk_numlines) + self.vtk_numline_components, self.vtk_numline_components):
+                for vtk_offsets in range(self.vtk_numline_components,
+                                        (self.vtk_numline_components * self.vtk_numlines) +
+                                        self.vtk_numline_components, self.vtk_numline_components):
                     f.write(pack('I', vtk_offsets))
                     pbar.update(n=4)
 
@@ -293,25 +303,24 @@ class GeometryView:
                 self.write_gprmax_info(f, G, materialsonly=True)
 
     def write_gprmax_info(self, f, G, materialsonly=False):
-        """
+        """Writes gprMax specific information relating material, source,
-        Writes gprMax specific information relating material, source,
+            and receiver names to numeric identifiers.
-        and receiver names to numeric identifiers.
 
         Args:
             f (filehandle): VTK file.
-            G (class): Grid class instance - holds essential parameters describing the model.
+            G (FDTDGrid): Parameters describing a grid in a model.
-            materialsonly (boolean): Only write information on materials
+            materialsonly (bool): Only write information on materials
         """
 
         f.write('\n\n<gprMax>\n'.encode('utf-8'))
         for material in G.materials:
-            f.write('<Material name="{}">{}</Material>\n'.format(material.ID, material.numID).encode('utf-8'))
+            f.write(f'<Material name="{material.ID}">{material.numID}</Material>\n'.encode('utf-8'))
         if not materialsonly:
             f.write('<PML name="PML boundary region">1</PML>\n'.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'))
+                f.write(f'<Sources name="{src.ID}">{index + 2}</Sources>\n'.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(f'<Receivers name="{rx.ID}">{index + 1}</Receivers>\n'.encode('utf-8'))
         f.write('</gprMax>\n'.encode('utf-8'))
 
 
@@ -334,8 +343,13 @@ class GeometryObjects:
         self.nx = self.xf - self.xs
         self.ny = self.yf - self.ys
         self.nz = self.zf - self.zs
-        self.filename = basefilename + '.h5'
+
-        self.materialsfilename = basefilename + '_materials.txt'
+        # Set filenames
+        parts = config.sim_config.input_file_path.with_suffix('').parts
+        self.filename_hdf5 = Path(*parts[:-1], basefilename)
+        self.filename_hdf5 = self.filename_hdf5.with_suffix('.h5')
+        self.filename_materials = Path(*parts[:-1], basefilename + '_materials')
+        self.filename_materials = self.filename_materials.with_suffix('.txt')
 
         # Sizes of arrays to write necessary to update progress bar
         self.solidsize = (self.nx + 1) * (self.ny + 1) * (self.nz + 1) * np.dtype(np.int16).itemsize
@@ -347,12 +361,12 @@ class GeometryObjects:
         """Write a geometry objects file in HDF5 format.
 
         Args:
-            G (class): Grid class instance - holds essential parameters describing the model.
+            G (FDTDGrid): Parameters describing a grid in a model.
             pbar (class): Progress bar class instance.
         """
 
         # Write the geometry objects to a HDF5 file
-        fdata = h5py.File(os.path.abspath(os.path.join(G.inputdirectory, self.filename)), 'w')
+        fdata = h5py.File(self.filename_hdf5, 'w')
         fdata.attrs['gprMax'] = __version__
         fdata.attrs['Title'] = G.title
         fdata.attrs['dx_dy_dz'] = (G.dx, G.dy, G.dz)
@@ -370,27 +384,28 @@ class GeometryObjects:
 
         # Write materials list to a text file
         # This includes all materials in range whether used in volume or not
-        fmaterials = open(os.path.abspath(os.path.join(G.inputdirectory, self.materialsfilename)), 'w')
+        fmaterials = open(self.filename_materials, 'w')
         for numID in range(minmat, maxmat + 1):
             for material in G.materials:
                 if material.numID == numID:
-                    fmaterials.write('#material: {:g} {:g} {:g} {:g} {}\n'.format(material.er, material.se, material.mr, material.sm, material.ID))
+                    fmaterials.write(f'#material: {material.er:g} {material.se:g} {material.mr:g} {material.sm:g} {material.ID}\n')
                     if material.poles > 0:
                         if 'debye' in material.type:
-                            dispersionstr = '#add_dispersion_debye: {:g} '.format(material.poles)
+                            dispersionstr = f'#add_dispersion_debye: {material.poles:g} '
                             for pole in range(material.poles):
-                                dispersionstr += '{:g} {:g} '.format(material.deltaer[pole], material.tau[pole])
+                                dispersionstr += f'{material.deltaer[pole]:g} {material.tau[pole]:g} '
                         elif 'lorenz' in material.type:
-                            dispersionstr = '#add_dispersion_lorenz: {:g} '.format(material.poles)
+                            dispersionstr = f'#add_dispersion_lorenz: {material.poles:g} '
                             for pole in range(material.poles):
-                                dispersionstr += '{:g} {:g} {:g} '.format(material.deltaer[pole], material.tau[pole], material.alpha[pole])
+                                dispersionstr += f'{material.deltaer[pole]:g} {material.tau[pole]:g} {material.alpha[pole]:g} '
                         elif 'drude' in material.type:
-                            dispersionstr = '#add_dispersion_drude: {:g} '.format(material.poles)
+                            dispersionstr = f'#add_dispersion_drude: {material.poles:g} '
                             for pole in range(material.poles):
-                                dispersionstr += '{:g} {:g} '.format(material.tau[pole], material.alpha[pole])
+                                dispersionstr += f'{material.tau[pole]:g} {material.alpha[pole]:g} '
                         dispersionstr += material.ID
                         fmaterials.write(dispersionstr + '\n')
 
+
 class GeometryViewFineMultiGrid:
     """Geometry view for all grids in the simulation.
 
@@ -400,11 +415,6 @@ class GeometryViewFineMultiGrid:
         view them at once in Paraview.
     """
 
-    if sys.byteorder == 'little':
-        byteorder = 'LittleEndian'
-    else:
-        byteorder = 'BigEndian'
-
     def __init__(self, xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, fileext, G):
         """
         Args:
@@ -432,7 +442,6 @@ class GeometryViewFineMultiGrid:
             # total additional lines required for subgrid
             self.additional_lines += sg_gv.n_total_lines
             # total additional points required for subgrid
-
             self.additional_points += sg_gv.n_total_points
 
         self.vtk_numpoints = self.additional_points + (self.nx + 1) * (self.ny + 1) * (self.nz + 1)
@@ -444,24 +453,21 @@ class GeometryViewFineMultiGrid:
         self.vtk_materials_offset = round_value(int(self.vtk_offsets_offset + (self.vtk_numlines * np.dtype(np.uint32).itemsize) + np.dtype(np.uint32).itemsize))
         self.datawritesize = np.dtype(np.float32).itemsize * self.vtk_numpoints * self.vtk_numpoint_components + np.dtype(np.uint32).itemsize * self.vtk_numlines * self.vtk_numline_components + np.dtype(np.uint32).itemsize * self.vtk_numlines + np.dtype(np.uint32).itemsize * self.vtk_numlines
 
-    def set_filename(self, appendmodelnumber):
+    def set_filename(self):
-        """
+        """Construct filename from user-supplied name and model run number."""
-        Construct filename from user-supplied name and model run number.
+        parts = config.model_configs[self.G.model_num].input_file_path.parts
-
+        self.filename = Path(*parts[:-1], parts[-1] + config.model_configs[self.G.model_num].appendmodelnumber)
-        Args:
+        self.filename = self.filename.with_suffix(self.fileext)
-            appendmodelnumber (str): Text to append to filename.
-        """
-
-        self.filename = os.path.abspath(os.path.join(os.path.dirname(os.path.abspath(config.general['inputfilepath'])), self.basefilename + appendmodelnumber))
-        self.filename += self.fileext
 
     def write_vtk(self, *args):
-        """
+        """Writes the geometry information to a VTK file.
-        Writes the geometry information to a VTK file. Either ImageData (.vti) for a
+            Either ImageData (.vti) for a per-cell geometry view, or
-        per-cell geometry view, or PolygonalData (.vtp) for a per-cell-edge geometry view.
+            PolygonalData (.vtp) for a per-cell-edge geometry view.
+
             N.B. No Python 3 support for VTK at time of writing (03/2015)
+
         Args:
-            G (class): Grid class instance - holds essential parameters describing the model.
+            G (FDTDGrid): Parameters describing a grid in a model.
         """
 
         G = self.G
@@ -471,7 +477,7 @@ class GeometryViewFineMultiGrid:
             # refine parameters for subgrid
 
             f.write('<?xml version="1.0"?>\n'.encode('utf-8'))
-            f.write('<VTKFile type="PolyData" version="1.0" byte_order="{}">\n'.format(GeometryViewFineMultiGrid.byteorder).encode('utf-8'))
+            f.write('<VTKFile type="PolyData" version="1.0" byte_order="{}">\n'.format(config.sim_config.vtk_byteorder).encode('utf-8'))
             f.write('<PolyData>\n<Piece NumberOfPoints="{}" NumberOfVerts="0" NumberOfLines="{}" NumberOfStrips="0" NumberOfPolys="0">\n'.format(self.vtk_numpoints, self.vtk_numlines).encode('utf-8'))
 
             f.write('<Points>\n<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="0" />\n</Points>\n'.encode('utf-8'))
@@ -484,7 +490,7 @@ class GeometryViewFineMultiGrid:
             f.write('</Piece>\n</PolyData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
 
             # Write points
-            print('writing points main grid')
+            log.info('writing points main grid')
             datasize = np.dtype(np.float32).itemsize * self.vtk_numpoints * self.vtk_numpoint_components
             f.write(pack('I', datasize))
             for i in range(0, G.nx + 1):
@@ -493,7 +499,7 @@ class GeometryViewFineMultiGrid:
                         f.write(pack('fff', i * G.dx, j * G.dy, k * G.dz))
 
             for sg_v in self.sg_views:
-                print('writing points subgrid')
+                log.info('writing points subgrid')
                 sg_v.write_points(f, G)
 
             n_x_lines = self.nx * (self.ny + 1) * (self.nz + 1)
@@ -508,7 +514,7 @@ class GeometryViewFineMultiGrid:
             z_lines = np.zeros((n_z_lines, 2), dtype=np.uint32)
             z_materials = np.zeros((n_z_lines), dtype=np.uint32)
 
-            print('calculate connectivity main grid')
+            log.info('calculate connectivity main grid')
             label = 0
             counter_x = 0
             counter_y = 0
@@ -540,7 +546,7 @@ class GeometryViewFineMultiGrid:
 
                         label = label + 1
 
-            print('calculate connectivity subgrids')
+            log.info('calculate connectivity subgrids')
             for sg_v in self.sg_views:
                 sg_v.populate_connectivity_and_materials(label)
                 # use the last subgrids label for the next view
@@ -584,13 +590,13 @@ class GeometryViewFineMultiGrid:
             #self.write_gprmax_info(f, G, materialsonly=True)
 
     def write_gprmax_info(self, f, G, materialsonly=False):
-        """
+        """Writes gprMax specific information relating material, source,
-        Writes gprMax specific information relating material, source,
+            and receiver names to numeric identifiers.
-        and receiver names to numeric identifiers.
+
         Args:
             f (filehandle): VTK file.
-            G (class): Grid class instance - holds essential parameters describing the model.
+            G (FDTDGrid): Parameters describing a grid in a model.
-            materialsonly (boolean): Only write information on materials
+            materialsonly (bool): Only write information on materials.
         """
 
         f.write('\n\n<gprMax>\n'.encode('utf-8'))
@@ -607,6 +613,7 @@ class GeometryViewFineMultiGrid:
 
 
 class SubgridGeometryView:
+
     def __init__(self, sg):
 
         self.sg = sg
@@ -646,8 +653,10 @@ class SubgridGeometryView:
                     f.write(pack('fff', p_x, p_y, p_z))
 
     def populate_connectivity_and_materials(self, label):
-        """Label is the starting label. 0 if no other grids are present but +1 the last label used
+        """Label is the starting label. 0 if no other grids are present but
-        for a multigrid view"""
+        +1 the last label used for a multigrid view.
+        """
+
         sg = self.sg
         self.label = label
 

gprMax/grid.py

@@ -174,21 +174,19 @@ class FDTDGrid:
         self.updatecoeffsH = np.zeros((len(self.materials), 5),
                                        dtype=config.sim_config.dtypes['float_or_double'])
 
-    def initialise_dispersive_arrays(self, dtype):
+    def initialise_dispersive_arrays(self):
-        """Initialise arrays for storing coefficients when there are dispersive materials present.
+        """Initialise arrays for storing coefficients when there are dispersive
-
+            materials present.
-        Args:
-            dtype (dtype): Dtype to use for dispersive arrays.
         """
         self.Tx = np.zeros((config.model_configs[self.model_num].materials['maxpoles'],
-                            self.nx + 1, self.ny + 1, self.nz + 1), dtype=dtype)
+                            self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.model_configs[self.model_num].materials['dispersivedtype'])
         self.Ty = np.zeros((config.model_configs[self.model_num].materials['maxpoles'],
-                            self.nx + 1, self.ny + 1, self.nz + 1), dtype=dtype)
+                            self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.model_configs[self.model_num].materials['dispersivedtype'])
         self.Tz = np.zeros((config.model_configs[self.model_num].materials['maxpoles'],
-                            self.nx + 1, self.ny + 1, self.nz + 1), dtype=dtype)
+                            self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.model_configs[self.model_num].materials['dispersivedtype'])
         self.updatecoeffsdispersive = np.zeros((len(self.materials), 3 *
                                                 config.model_configs[self.model_num].materials['maxpoles']),
-                                                dtype=dtype)
+                                                dtype=config.model_configs[self.model_num].materials['dispersivedtype'])
 
     def mem_est_basic(self):
         """Estimate the amount of memory (RAM) required for grid arrays.
@@ -275,6 +273,8 @@ class FDTDGrid:
         """Clear arrays for field components and PMLs."""
         # Clear arrays for field components
         self.initialise_field_arrays()
+        if config.model_configs[self.model_num].materials['maxpoles'] != 0:
+            self.initialise_dispersive_arrays()
 
         # Clear arrays for fields in PML
         for pml in self.pmls:
@@ -313,12 +313,19 @@ class CUDAGrid(FDTDGrid):
         self.bpg = (int(np.ceil(((self.nx + 1) * (self.ny + 1) *
                    (self.nz + 1)) / self.tpb[0])), 1, 1)
 
-    def initialise_arrays(self):
+    def initialise_geometry_arrays(self):
+        """Initialise an array for cell edge IDs (ID) on GPU."""
+        import pycuda.gpuarray as gpuarray
+
+        super().initialise_geometry_arrays()
+        self.ID_gpu = gpuarray.to_gpu(self.ID)
+
+    def initialise_field_arrays(self):
         """Initialise geometry and field arrays on GPU."""
 
         import pycuda.gpuarray as gpuarray
 
-        self.ID_gpu = gpuarray.to_gpu(self.ID)
+        super().initialise_field_arrays()
         self.Ex_gpu = gpuarray.to_gpu(self.Ex)
         self.Ey_gpu = gpuarray.to_gpu(self.Ey)
         self.Ez_gpu = gpuarray.to_gpu(self.Ez)
@@ -326,22 +333,31 @@ class CUDAGrid(FDTDGrid):
         self.Hy_gpu = gpuarray.to_gpu(self.Hy)
         self.Hz_gpu = gpuarray.to_gpu(self.Hz)
 
-    def initialise_dispersive_arrays(self, dtype):
+    def initialise_dispersive_arrays(self):
-        """Initialise dispersive material coefficient arrays on GPU.
+        """Initialise dispersive material coefficient arrays on GPU."""
-
-        Args:
-            dtype (dtype): Dtype to use for dispersive arrays.
-        """
 
         import pycuda.gpuarray as gpuarray
 
-        super().initialise_dispersive_arrays(dtype)
+        super().initialise_dispersive_arrays()
-
         self.Tx_gpu = gpuarray.to_gpu(self.Tx)
         self.Ty_gpu = gpuarray.to_gpu(self.Ty)
         self.Tz_gpu = gpuarray.to_gpu(self.Tz)
         self.updatecoeffsdispersive_gpu = gpuarray.to_gpu(self.updatecoeffsdispersive)
 
+    def reset_fields(self):
+        """Clear arrays for field components and PMLs."""
+        
+        super().reset_fields()
+
+        # Clear arrays for field components
+        self.initialise_field_arrays()
+        if config.model_configs[self.model_num].materials['maxpoles'] != 0:
+            self.initialise_dispersive_arrays()
+
+        # Clear arrays for fields in PML
+        for pml in self.pmls:
+            pml.initialise_field_arrays()
+
     def memory_check(self, snapsmemsize=0):
         """Check if model can be run on specified GPU."""
 
@@ -487,7 +503,7 @@ def Ix(x, y, z, Hx, Hy, Hz, G):
     Args:
         x, y, z (float): Coordinates of position in grid.
         Hx, Hy, Hz (memory view): numpy array of magnetic field values.
-        G (FDTDGrid): Holds essential parameters describing a model.
+        G (FDTDGrid): Parameters describing a grid in a model.
     """
 
     if y == 0 or z == 0:
@@ -504,7 +520,7 @@ def Iy(x, y, z, Hx, Hy, Hz, G):
     Args:
         x, y, z (float): Coordinates of position in grid.
         Hx, Hy, Hz (memory view): numpy array of magnetic field values.
-        G (FDTDGrid): Holds essential parameters describing a model.
+        G (FDTDGrid): Parameters describing a grid in a model.
     """
 
     if x == 0 or z == 0:
@@ -521,7 +537,7 @@ def Iz(x, y, z, Hx, Hy, Hz, G):
     Args:
         x, y, z (float): Coordinates of position in grid.
         Hx, Hy, Hz (memory view): numpy array of magnetic field values.
-        G (FDTDGrid): Holds essential parameters describing a model.
+        G (FDTDGrid): Parameters describing a grid in a model.
     """
 
     if x == 0 or y == 0:

gprMax/hash_cmds_file.py

@@ -19,6 +19,7 @@
 from io import StringIO
 import logging
 import os
+from pathlib import Path
 import sys
 
 import gprMax.config as config
@@ -157,23 +158,22 @@ def process_include_files(hashcmds, inputfile):
     return processedincludecmds
 
 
-def write_processed_file(processedlines, appendmodelnumber):
+def write_processed_file(processedlines, G):
     """Writes an input file after any Python code and include commands
         in the original input file have been processed.
 
     Args:
         processedlines (list): Input commands after after processing any
             Python code and include commands.
-        appendmodelnumber (str): Text to append to filename.
+        G (FDTDGrid): Parameters describing a grid in a model.
     """
 
-    processedfile = (os.path.join(config.general['outputfilepath'],
+    parts = config.model_configs[G.model_num].output_file_path.parts
-                     os.path.splitext(config.general['inputfilepath'])[0] +
+    processedfile = (Path(*parts[:-1], parts[-1] + '_processed.in'))
-                     appendmodelnumber + '_processed.in'))
 
     with open(processedfile, 'w') as f:
         for item in processedlines:
-            f.write('{}'.format(item))
+            f.write(f'{item}')
 
     log.info(f'Written input commands, after processing any Python code and include commands, to file: {processedfile}\n')
 
@@ -282,13 +282,12 @@ def get_user_objects(processedlines, check=True):
     return user_objs
 
 
-def parse_hash_commands(model_config, G, scene):
+def parse_hash_commands(scene, G):
     """Parse user hash commands and add them to the scene.
 
     Args:
-        model_config (ModelConfig): Model level configuration object.
-        G (FDTDGrid): Holds essential parameters describing a model.
         scene (Scene): Scene object.
+        G (FDTDGrid): Parameters describing a grid in a model.
 
     Returns:
         scene (Scene): Scene object.
@@ -296,7 +295,7 @@ def parse_hash_commands(model_config, G, scene):
 
     with open(config.sim_config.input_file_path) as inputfile:
 
-        usernamespace = model_config.get_usernamespace()
+        usernamespace = config.model_configs[G.model_num].get_usernamespace()
 
         # Read input file and process any Python and include file commands
         processedlines = process_python_include_code(inputfile, usernamespace)
@@ -311,7 +310,7 @@ def parse_hash_commands(model_config, G, scene):
         # Write a file containing the input commands after Python or include
         # file commands have been processed
         if config.sim_config.args.write_processed:
-            write_processed_file(processedlines, model_config.appendmodelnumber, G)
+            write_processed_file(processedlines, G)
 
         user_objs = get_user_objects(processedlines, check=True)
         for user_obj in user_objs:

gprMax/materials.py

@@ -49,7 +49,7 @@ class Material:
         """Calculates the magnetic update coefficients of the material.
 
         Args:
-            G (FDTDGrid): Holds essential parameters describing a model.
+            G (FDTDGrid): Parameters describing a grid in a model.
         """
 
         HA = (config.m0 * self.mr / G.dt) + 0.5 * self.sm
@@ -64,7 +64,7 @@ class Material:
         """Calculates the electric update coefficients of the material.
 
         Args:
-            G (FDTDGrid): Holds essential parameters describing a model.
+            G (FDTDGrid): Parameters describing a grid in a model.
         """
 
         EA = (config.sim_config.em_consts['e0'] * self.er / G.dt) + 0.5 * self.se
@@ -126,7 +126,7 @@ class DispersiveMaterial(Material):
         """Calculates the electric update coefficients of the material.
 
         Args:
-            G (FDTDGrid): Holds essential parameters describing a model.
+            G (FDTDGrid): Parameters describing a grid in a model.
         """
 
         # The implementation of the dispersive material modelling comes from the
@@ -216,7 +216,7 @@ def process_materials(G):
         store in arrays, and build text list of materials/properties
 
     Args:
-        G (FDTDGrid): Holds essential parameters describing a model.
+        G (FDTDGrid): Parameters describing a grid in a model.
 
     Returns:
         materialsdata (list): List of material IDs, names, and properties to
@@ -315,7 +315,7 @@ class PeplinskiSoil:
 
         Args:
             nbins (int): Number of bins to use to create the different materials.
-            G (FDTDGrid): Holds essential parameters describing a model.
+            G (FDTDGrid): Parameters describing a grid in a model.
         """
 
         # Debye model properties of water
@@ -382,7 +382,7 @@ def create_built_in_materials(G):
     """Create pre-defined (built-in) materials.
 
     Args:
-        G (FDTDGrid): Holds essential parameters describing a model.
+        G (FDTDGrid): Parameters describing a grid in a model.
     """
 
     G.n_built_in_materials = len(G.materials)

gprMax/model_build_run.py

@@ -103,6 +103,8 @@ class ModelBuildRun:
         # Normal model reading/building process; bypassed if geometry information to be reused
         self.reuse_geometry() if config.model_configs[G.model_num].reuse_geometry else self.build_geometry()
 
+        log.info(f'\nOutput path: {config.model_configs[G.model_num].output_file_path.parent}')
+
         # Adjust position of simple sources and receivers if required
         if G.srcsteps[0] != 0 or G.srcsteps[1] != 0 or G.srcsteps[2] != 0:
             for source in itertools.chain(G.hertziandipoles, G.magneticdipoles):
@@ -135,7 +137,8 @@ class ModelBuildRun:
         if not (G.geometryviews or G.geometryobjectswrite) and config.sim_config.args.geometry_only:
             log.warning(Fore.RED + f'\nNo geometry views or geometry objects found.' + Style.RESET_ALL)
         for i, geometryview in enumerate(G.geometryviews):
-            geometryview.set_filename(config.model_configs[G.model_num].appendmodelnumber)
+            log.info('')
+            geometryview.set_filename()
             pbar = tqdm(total=geometryview.datawritesize, unit='byte', unit_scale=True,
                         desc=f'Writing geometry view file {i + 1}/{len(G.geometryviews)}, {geometryview.filename.name}',
                         ncols=get_terminal_width() - 1, file=sys.stdout,
@@ -143,6 +146,7 @@ class ModelBuildRun:
             geometryview.write_vtk(G, pbar)
             pbar.close()
         for i, geometryobject in enumerate(G.geometryobjectswrite):
+            log.info('')
             pbar = tqdm(total=geometryobject.datawritesize, unit='byte', unit_scale=True,
                         desc=f'Writing geometry object file {i + 1}/{len(G.geometryobjectswrite)}, {geometryobject.filename.name}',
                         ncols=get_terminal_width() - 1, file=sys.stdout,
@@ -157,22 +161,11 @@ class ModelBuildRun:
 
         scene = self.build_scene()
 
-        # Combine available grids and check memory requirements
+        # Combine available grids and check basic memory requirements
         grids = [G] + G.subgrids
         for grid in grids:
             config.model_configs[G.model_num].mem_use += grid.mem_est_basic()
         mem_check(config.model_configs[G.model_num].mem_use)
-        log.info(f'\nMemory (RAM) required: ~{human_size(config.model_configs[G.model_num].mem_use)}')
-
-        gridbuilders = [GridBuilder(grid) for grid in grids]
-
-        for gb in gridbuilders:
-            pml_information(gb.grid)
-            gb.build_pmls()
-            gb.build_components()
-            gb.tm_grid_update()
-            gb.update_voltage_source_materials()
-            gb.grid.initialise_std_update_coeff_arrays()
 
         # Set datatype for dispersive arrays if there are any dispersive materials.
         if config.model_configs[G.model_num].materials['maxpoles'] != 0:
@@ -187,10 +180,6 @@ class ModelBuildRun:
             # Update estimated memory (RAM) usage
             config.model_configs[G.model_num].mem_use += G.mem_est_dispersive()
             mem_check(config.model_configs[G.model_num].mem_use)
-            log.info(f'Memory (RAM) required - updated (dispersive): ~{human_size(config.model_configs[G.model_num].mem_use)}')
-
-            for gb in gridbuilders:
-                gb.grid.initialise_dispersive_arrays(config.model_configs[G.model_num].materials['dispersivedtype'])
 
         # Check there is sufficient memory to store any snapshots
         if G.snapshots:
@@ -203,10 +192,20 @@ class ModelBuildRun:
             mem_check(config.model_configs[G.model_num].mem_use)
             if config.sim_config.general['cuda']:
                 mem_check_gpu_snaps(G.model_num, snaps_mem)
-            log.info(f'Memory (RAM) required - updated (snapshots): ~{human_size(config.model_configs[G.model_num].mem_use)}')
 
-        # Build materials
+        log.info(f'\nMemory (RAM) required: ~{human_size(config.model_configs[G.model_num].mem_use)}')
+
+        # Build grids
+        gridbuilders = [GridBuilder(grid) for grid in grids]
         for gb in gridbuilders:
+            pml_information(gb.grid)
+            gb.build_pmls()
+            gb.build_components()
+            gb.tm_grid_update()
+            gb.update_voltage_source_materials()
+            gb.grid.initialise_std_update_coeff_arrays()
+            if config.model_configs[G.model_num].materials['maxpoles'] != 0:
+                gb.grid.initialise_dispersive_arrays()
             gb.build_materials()
 
         # Check to see if numerical dispersion might be a problem
@@ -237,25 +236,13 @@ class ModelBuildRun:
         if not scene:
             scene = Scene()
             # Parse the input file into user objects and add them to the scene
-            scene = parse_hash_commands(config.model_configs[self.G.model_num], self.G, scene)
+            scene = parse_hash_commands(scene, self.G)
 
         # Creates the internal simulation objects
         scene.create_internal_objects(self.G)
 
         return scene
 
-    def create_output_directory(self):
-        log.debug('Fix output directory path setting')
-        # if self.G.outputdirectory:
-        #     # Check and set output directory and filename
-        #     try:
-        #         os.mkdir(self.G.outputdirectory)
-        #         log.info(f'\nCreated output directory: {self.G.outputdirectory}')
-        #     except FileExistsError:
-        #         pass
-        #     # Modify the output path (hack)
-        #     config.model_configs[G.model_num].output_file_path_ext = Path(self.G.outputdirectory, config.model_configs[G.model_num].output_file_path_ext)
-
     def write_output_data(self):
         """Write output data, i.e. field data for receivers and snapshots
             to file(s).
@@ -266,19 +253,17 @@ class ModelBuildRun:
 
         # Write any snapshots to file
         if self.G.snapshots:
-            # Create directory and construct filename from user-supplied name
+            # Create directory for snapshots
-            # and model run number
+            config.model_configs[self.G.model_num].set_snapshots_file_path()
-            snapshotdir = config.model_configs[self.G.model_num].snapshot_dir
+            snapshotdir = config.model_configs[self.G.model_num].snapshot_file_path
-            if not os.path.exists(snapshotdir):
+            snapshotdir.mkdir(exist_ok=True)
-                os.mkdir(snapshotdir)
 
             log.info('')
             for i, snap in enumerate(self.G.snapshots):
-                fn = snapshotdir / Path(config.model_configs[self.G.model_num].output_file_path.stem + '_' + snap.basefilename)
+                fn = snapshotdir / Path(snap.filename)
                 snap.filename = fn.with_suffix('.vti')
                 pbar = tqdm(total=snap.vtkdatawritesize, leave=True, unit='byte',
-                            unit_scale=True, desc=f'Writing snapshot file {i + 1} of {len(self.G.snapshots)}, {os.path.split(snap.filename)[1]}', ncols=get_terminal_width() - 1, file=sys.stdout,
+                            unit_scale=True, desc=f'Writing snapshot file {i + 1} of {len(self.G.snapshots)}, {snap.filename.name}', ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.sim_config.general['progressbars'])
-                            disable=not config.general['progressbars'])
                 snap.write_vtk_imagedata(pbar, self.G)
                 pbar.close()
             log.info('')
@@ -308,7 +293,6 @@ class ModelBuildRun:
             tsolve (float): time taken to execute solving (seconds).
         """
 
-        self.create_output_directory()
         log.info(f'\nOutput file: {config.model_configs[self.G.model_num].output_file_path_ext}')
 
         # Check number of OpenMP threads

gprMax/pml.py

@@ -339,11 +339,13 @@ class CUDAPML(PML):
         solving on GPU using CUDA.
     """
 
-    def initialise_field_arrays_gpu(self):
+    def initialise_field_arrays(self):
         """Initialise PML field and coefficient arrays on GPU."""
 
         import pycuda.gpuarray as gpuarray
 
+        super().initialise_field_arrays()
+
         self.ERA_gpu = gpuarray.to_gpu(self.ERA)
         self.ERB_gpu = gpuarray.to_gpu(self.ERB)
         self.ERE_gpu = gpuarray.to_gpu(self.ERE)
@@ -352,52 +354,16 @@ class CUDAPML(PML):
         self.HRB_gpu = gpuarray.to_gpu(self.HRB)
         self.HRE_gpu = gpuarray.to_gpu(self.HRE)
         self.HRF_gpu = gpuarray.to_gpu(self.HRF)
-
+        self.EPhi1_gpu = gpuarray.to_gpu(self.EPhi1)
-        if self.direction[0] == 'x':
+        self.EPhi2_gpu = gpuarray.to_gpu(self.EPhi2)
-            self.EPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS),
+        self.HPhi1_gpu = gpuarray.to_gpu(self.HPhi1)
-                                             self.nx + 1, self.ny, self.nz + 1),
+        self.HPhi2_gpu = gpuarray.to_gpu(self.HPhi2)
-                                             dtype=config.sim_config.dtypes['float_or_double']))
-            self.EPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS),
-                                             self.nx + 1, self.ny + 1, self.nz),
-                                             dtype=config.sim_config.dtypes['float_or_double']))
-            self.HPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS),
-                                             self.nx, self.ny + 1, self.nz),
-                                             dtype=config.sim_config.dtypes['float_or_double']))
-            self.HPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS),
-                                             self.nx, self.ny, self.nz + 1),
-                                             dtype=config.sim_config.dtypes['float_or_double']))
-        elif self.direction[0] == 'y':
-            self.EPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS),
-                                             self.nx, self.ny + 1, self.nz + 1),
-                                             dtype=config.sim_config.dtypes['float_or_double']))
-            self.EPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS),
-                                             self.nx + 1, self.ny + 1, self.nz),
-                                             dtype=config.sim_config.dtypes['float_or_double']))
-            self.HPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS),
-                                             self.nx + 1, self.ny, self.nz),
-                                             dtype=config.sim_config.dtypes['float_or_double']))
-            self.HPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS),
-                                             self.nx, self.ny, self.nz + 1),
-                                             dtype=config.sim_config.dtypes['float_or_double']))
-        elif self.direction[0] == 'z':
-            self.EPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS),
-                                             self.nx, self.ny + 1, self.nz + 1),
-                                             dtype=config.sim_config.dtypes['float_or_double']))
-            self.EPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS),
-                                             self.nx + 1, self.ny, self.nz + 1),
-                                             dtype=config.sim_config.dtypes['float_or_double']))
-            self.HPhi1_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS),
-                                             self.nx + 1, self.ny, self.nz),
-                                             dtype=config.sim_config.dtypes['float_or_double']))
-            self.HPhi2_gpu = gpuarray.to_gpu(np.zeros((len(self.CFS),
-                                             self.nx, self.ny + 1, self.nz),
-                                             dtype=config.sim_config.dtypes['float_or_double']))
 
     def set_blocks_per_grid(self, G):
         """Set the blocks per grid size used for updating the PML field arrays on a GPU.
 
         Args:
-            G (FDTDGrid): Holds essential parameters describing the model.
+            G (FDTDGrid): Parameters describing a grid in a model.
         """
 
         self.bpg = (int(np.ceil(((self.EPhi1_gpu.shape[1] + 1) *
@@ -422,7 +388,7 @@ class CUDAPML(PML):
             correction on the GPU.
 
         Args:
-            G (FDTDGrid): Holds essential parameters describing the model.
+            G (FDTDGrid): Parameters describing a grid in a model.
         """
 
         self.update_electric_gpu(np.int32(self.xs), np.int32(self.xf),
@@ -449,7 +415,7 @@ class CUDAPML(PML):
             correction on the GPU.
 
         Args:
-            G (FDTDGrid): Holds essential parameters describing the model.
+            G (FDTDGrid): Parameters describing a grid in a model.
         """
         self.update_magnetic_gpu(np.int32(self.xs), np.int32(self.xf),
                                  np.int32(self.ys), np.int32(self.yf),
@@ -474,7 +440,7 @@ def pml_information(G):
     """Information about PMLs.
 
     Args:
-        G (FDTDGrid): Holds essential parameters describing the model.
+        G (FDTDGrid): Parameters describing a grid in a model.
     """
     # No PML
     if all(value == 0 for value in G.pmlthickness.values()):
@@ -497,7 +463,7 @@ def build_pml(G, key, value):
         (based on underlying material er and mr from solid array).
 
     Args:
-        G (FDTDGrid): Holds essential parameters describing the model.
+        G (FDTDGrid): Parameters describing a grid in a model.
         key (str): Identifier of PML slab.
         value (int): Thickness of PML slab in cells.
     """

gprMax/snapshots.py

@@ -39,11 +39,6 @@ class Snapshot:
     # GPU - blocks per grid - set according to largest requested snapshot
     bpg = None
 
-    # Set string for byte order
-    byteorder = 'LittleEndian' if sys.byteorder == 'little' else 'BigEndian'
-
-
-
     def __init__(self, xs=None, ys=None, zs=None, xf=None, yf=None, zf=None,
                  dx=None, dy=None, dz=None, time=None, filename=None):
         """
@@ -54,14 +49,6 @@ class Snapshot:
             filename (str): Filename to save to.
         """
 
-        # Set format text and string depending on float type
-        if config.sim_config.dtypes['float_or_double'] == np.float32:
-            self.floatname = 'Float32'
-            self.floatstring = 'f'
-        elif config.sim_config.dtypes['float_or_double'] == np.float64:
-            self.floatname = 'Float64'
-            self.floatstring = 'd'
-
         self.fieldoutputs = {'electric': True, 'magnetic': True}
         self.xs = xs
         self.ys = ys
@@ -79,21 +66,20 @@ class Snapshot:
         self.sy = slice(self.ys, self.yf + self.dy, self.dy)
         self.sz = slice(self.zs, self.zf + self.dz, self.dz)
         self.ncells = self.nx * self.ny * self.nz
-        self.datasizefield = (3 * np.dtype(config.dtypes['float_or_double']).itemsize
+        self.datasizefield = (3 * np.dtype(config.sim_config.dtypes['float_or_double']).itemsize
                                 * self.ncells)
-        self.vtkdatawritesize = ((self.fieldoutputs['electric']
+        self.vtkdatawritesize = ((self.fieldoutputs['electric'] +
-                                  + self.fieldoutputs['magnetic']) * self.datasizefield
+                                  self.fieldoutputs['magnetic']) *
-                                 + (self.fieldoutputs['electric']
+                                  self.datasizefield + (self.fieldoutputs['electric'] +
-                                    + self.fieldoutputs['magnetic'])
+                                  self.fieldoutputs['magnetic']) * np.dtype(np.uint32).itemsize)
-                                 * np.dtype(np.uint32).itemsize)
         self.time = time
-        self.basefilename = filename
+        self.filename = filename
 
     def store(self, G):
         """Store (in memory) electric and magnetic field values for snapshot.
 
         Args:
-            G (class): Grid class instance - holds essential parameters describing the model.
+            G (FDTDGrid): Parameters describing a grid in a model.
         """
 
         # Memory views of field arrays to dimensions required for the snapshot
@@ -105,12 +91,12 @@ class Snapshot:
         Hzslice = np.ascontiguousarray(G.Hz[self.sx, self.sy, self.sz])
 
         # Create arrays to hold the field data for snapshot
-        Exsnap = np.zeros((self.nx, self.ny, self.nz), dtype=config.dtypes['float_or_double'])
+        Exsnap = np.zeros((self.nx, self.ny, self.nz), dtype=config.sim_config.dtypes['float_or_double'])
-        Eysnap = np.zeros((self.nx, self.ny, self.nz), dtype=config.dtypes['float_or_double'])
+        Eysnap = np.zeros((self.nx, self.ny, self.nz), dtype=config.sim_config.dtypes['float_or_double'])
-        Ezsnap = np.zeros((self.nx, self.ny, self.nz), dtype=config.dtypes['float_or_double'])
+        Ezsnap = np.zeros((self.nx, self.ny, self.nz), dtype=config.sim_config.dtypes['float_or_double'])
-        Hxsnap = np.zeros((self.nx, self.ny, self.nz), dtype=config.dtypes['float_or_double'])
+        Hxsnap = np.zeros((self.nx, self.ny, self.nz), dtype=config.sim_config.dtypes['float_or_double'])
-        Hysnap = np.zeros((self.nx, self.ny, self.nz), dtype=config.dtypes['float_or_double'])
+        Hysnap = np.zeros((self.nx, self.ny, self.nz), dtype=config.sim_config.dtypes['float_or_double'])
-        Hzsnap = np.zeros((self.nx, self.ny, self.nz), dtype=config.dtypes['float_or_double'])
+        Hzsnap = np.zeros((self.nx, self.ny, self.nz), dtype=config.sim_config.dtypes['float_or_double'])
 
         # Calculate field values at points (comes from averaging field components in cells)
         calculate_snapshot_fields(
@@ -141,28 +127,28 @@ class Snapshot:
 
         Args:
             pbar (class): Progress bar class instance.
-            G (class): Grid class instance - holds essential parameters describing the model.
+            G (FDTDGrid): Parameters describing a grid in a model.
         """
 
-        hfield_offset = (3 * np.dtype(config.dtypes['float_or_double']).itemsize
+        hfield_offset = (3 * np.dtype(config.sim_config.dtypes['float_or_double']).itemsize
                          * self.ncells + np.dtype(np.uint32).itemsize)
 
         self.filehandle = open(self.filename, 'wb')
         self.filehandle.write('<?xml version="1.0"?>\n'.encode('utf-8'))
-        self.filehandle.write(f'<VTKFile type="ImageData" version="1.0" byte_order="{Snapshot.byteorder}">\n'.encode('utf-8'))
+        self.filehandle.write(f'<VTKFile type="ImageData" version="1.0" byte_order="{config.sim_config.vtk_byteorder}">\n'.encode('utf-8'))
         self.filehandle.write(f'<ImageData WholeExtent="{self.xs} {round_value(self.xf / self.dx)} {self.ys} {round_value(self.yf / self.dy)} {self.zs} {round_value(self.zf / self.dz)}" Origin="0 0 0" Spacing="{self.dx * G.dx:.3} {self.dy * G.dy:.3} {self.dz * G.dz:.3}">\n'.encode('utf-8'))
         self.filehandle.write(f'<Piece Extent="{self.xs} {round_value(self.xf / self.dx)} {self.ys} {round_value(self.yf / self.dy)} {self.zs} {round_value(self.zf / self.dz)}">\n'.encode('utf-8'))
 
         if self.fieldoutputs['electric'] and self.fieldoutputs['magnetic']:
             self.filehandle.write('<CellData Vectors="E-field H-field">\n'.encode('utf-8'))
-            self.filehandle.write(f'<DataArray type="{self.floatname}" Name="E-field" NumberOfComponents="3" format="appended" offset="0" />\n'.encode('utf-8'))
+            self.filehandle.write(f"""<DataArray type="{config.sim_config.dtypes['vtk_float']}" Name="E-field" NumberOfComponents="3" format="appended" offset="0" />\n""".encode('utf-8'))
-            self.filehandle.write(f'<DataArray type="{self.floatname}" Name="H-field" NumberOfComponents="3" format="appended" offset="{hfield_offset}" />\n'.encode('utf-8'))
+            self.filehandle.write(f"""<DataArray type="{config.sim_config.dtypes['vtk_float']}" Name="H-field" NumberOfComponents="3" format="appended" offset="{hfield_offset}" />\n""".encode('utf-8'))
         elif self.fieldoutputs['electric']:
             self.filehandle.write('<CellData Vectors="E-field">\n'.encode('utf-8'))
-            self.filehandle.write(f'<DataArray type="{self.floatname}" Name="E-field" NumberOfComponents="3" format="appended" offset="0" />\n'.encode('utf-8'))
+            self.filehandle.write(f"""<DataArray type="{config.sim_config.dtypes['vtk_float']}" Name="E-field" NumberOfComponents="3" format="appended" offset="0" />\n""".encode('utf-8'))
         elif self.fieldoutputs['magnetic']:
             self.filehandle.write('<CellData Vectors="H-field">\n'.encode('utf-8'))
-            self.filehandle.write(f'<DataArray type="{self.floatname}" Name="H-field" NumberOfComponents="3" format="appended" offset="0" />\n'.encode('utf-8'))
+            self.filehandle.write(f"""<DataArray type="{config.sim_config.dtypes['vtk_float']}" Name="H-field" NumberOfComponents="3" format="appended" offset="0" />\n""".encode('utf-8'))
 
         self.filehandle.write('</CellData>\n</Piece>\n</ImageData>\n<AppendedData encoding="raw">\n_'.encode('utf-8'))
 
@@ -188,7 +174,7 @@ def initialise_snapshot_array_gpu(G):
     """Initialise array on GPU for to store field data for snapshots.
 
     Args:
-        G (class): Grid class instance - holds essential parameters describing the model.
+        G (FDTDGrid): Parameters describing a grid in a model.
 
     Returns:
         snapE_gpu, snapH_gpu (float): numpy arrays of snapshot data on GPU.
@@ -206,24 +192,26 @@ def initialise_snapshot_array_gpu(G):
             Snapshot.nz_max = snap.nz
 
     # GPU - blocks per grid - according to largest requested snapshot
-    Snapshot.bpg = (int(np.ceil(((Snapshot.nx_max) * (Snapshot.ny_max) * (Snapshot.nz_max)) / Snapshot.tpb[0])), 1, 1)
+    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);
     # if snapshots are not being stored on the GPU during the simulation then
     # they are copied back to the host after each iteration, hence numsnaps = 1
-    numsnaps = 1 if config.cuda['snapsgpu2cpu'] else len(G.snapshots)
+    numsnaps = 1 if config.model_configs[G.model_num].cuda['snapsgpu2cpu'] else len(G.snapshots)
     snapEx = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
-                      dtype=config.dtypes['float_or_double'])
+                      dtype=config.sim_config.dtypes['float_or_double'])
     snapEy = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
-                      dtype=config.dtypes['float_or_double'])
+                      dtype=config.sim_config.dtypes['float_or_double'])
     snapEz = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
-                      dtype=config.dtypes['float_or_double'])
+                      dtype=config.sim_config.dtypes['float_or_double'])
     snapHx = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
-                      dtype=config.dtypes['float_or_double'])
+                      dtype=config.sim_config.dtypes['float_or_double'])
     snapHy = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
-                      dtype=config.dtypes['float_or_double'])
+                      dtype=config.sim_config.dtypes['float_or_double'])
     snapHz = np.zeros((numsnaps, Snapshot.nx_max, Snapshot.ny_max, Snapshot.nz_max),
-                      dtype=config.dtypes['float_or_double'])
+                      dtype=config.sim_config.dtypes['float_or_double'])
 
     # Copy arrays to GPU
     snapEx_gpu = gpuarray.to_gpu(snapEx)

gprMax/subgrids/updates.py

@@ -75,7 +75,7 @@ class SubgridUpdater(CPUUpdates):
             subgrid (SubGrid3d): Subgrid to be updated.
             precursors (PrecursorNodes): Precursor nodes associated with
             the subgrid - contain interpolated fields.
-            G (FDTDGrid): Holds essential parameters describing a model.
+            G (FDTDGrid): Parameters describing a grid in a model.
         """
         super().__init__(subgrid)
         self.precursors = precursors