diff --git a/gprMax/cmds_geometry/geometry_objects_read.py b/gprMax/cmds_geometry/geometry_objects_read.py
index c00d7386..8bfafd36 100644
--- a/gprMax/cmds_geometry/geometry_objects_read.py
+++ b/gprMax/cmds_geometry/geometry_objects_read.py
@@ -93,7 +93,7 @@ class GeometryObjectsRead(UserObjectGeometry):
         f = h5py.File(geofile, 'r')
         dx_dy_dz = f.attrs['dx_dy_dz']
         if round_value(dx_dy_dz[0] / G.dx) != 1 or round_value(dx_dy_dz[1] / G.dy) != 1 or round_value(dx_dy_dz[2] / G.dz) != 1:
-            raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires the spatial resolution of the geometry objects file to match the spatial resolution of the model')
+            raise CmdInputError(self.__str__() + ' requires the spatial resolution of the geometry objects file to match the spatial resolution of the model')
 
         data = f['/data'][:]
 
diff --git a/gprMax/config.py b/gprMax/config.py
index d890c664..dbc8f684 100644
--- a/gprMax/config.py
+++ b/gprMax/config.py
@@ -40,7 +40,7 @@ log = logging.getLogger(__name__)
 # Single instance of SimConfig to hold simulation configuration parameters.
 sim_config = None
 
-# Instance of ModelConfig that hold model configuration parameters.
+# Instances of ModelConfig that hold model configuration parameters.
 model_configs = []
 
 # Each model in a simulation is given a unique number when the instance of
@@ -74,24 +74,21 @@ class ModelConfig:
                          'snapsgpu2cpu': False}
 
         # Total memory usage for all grids in the model. Starts with 50MB overhead.
-        self.mem_use = 50e6
+        self.mem_overhead = 50e6
+        self.mem_use = self.mem_overhead
 
         self.reuse_geometry = False
 
         # String to print at start of each model run
-        inputfilestr = f'\n--- Model {model_num + 1}/{sim_config.model_end}, input file: {sim_config.input_file_path}'
-        self.set_inputfilestr(inputfilestr)
+        s = f'\n--- Model {model_num + 1}/{sim_config.model_end}, input file: {sim_config.input_file_path}'
+        self.inputfilestr = Fore.GREEN + f"{s} {'-' * (get_terminal_width() - 1 - len(s))}\n" + Style.RESET_ALL
 
-        if not sim_config.single_model:
-            self.appendmodelnumber = str(model_num + 1) # Indexed from 1
-        else:
-            self.appendmodelnumber = ''
-
-        # Output file path for specific model
+        # Output file path and name for specific model
+        self.appendmodelnumber = '' if sim_config.single_model else str(model_num + 1) # Indexed from 1
         self.set_output_file_path()
 
         # Specify a snapshot directory
-        self.set_snapshots_file_path()
+        self.set_snapshots_dir()
 
         # Numerical dispersion analysis parameters
         #   highestfreqthres: threshold (dB) down from maximum power (0dB) of main frequency used
@@ -124,13 +121,18 @@ class ModelConfig:
                 'current_model_run': model_num + 1,
                 'inputfile': sim_config.input_file_path.resolve()}
 
-    def set_inputfilestr(self, inputfilestr):
-        """Set string describing model.
-
-        Args:
-            inputfilestr (str): Description of model.
+    def set_dispersive_material_types(self):
+        """Set data type for disperive materials. Complex if Drude or Lorentz
+            materials are present. Real if Debye materials.
         """
-        self.inputfilestr = Fore.GREEN + f"{inputfilestr} {'-' * (get_terminal_width() - 1 - len(inputfilestr))}\n" + Style.RESET_ALL
+        if self.materials['drudelorentz']:
+            self.materials['cuda_real_func'] = '.real()'
+            self.materials['dispersivedtype'] = sim_config.dtypes['complex']
+            self.materials['dispersiveCdtype'] = sim_config.dtypes['C_complex']
+        else:
+            self.materials['cuda_real_func'] = ''
+            self.materials['dispersivedtype'] = sim_config.dtypes['float_or_double']
+            self.materials['dispersiveCdtype'] = sim_config.dtypes['C_float_or_double']
 
     def set_output_file_path(self, outputdir=None):
         """Output file path can be provided by the user via the API or an input file
@@ -156,10 +158,11 @@ class ModelConfig:
         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):
+    def set_snapshots_dir(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')
+        self.snapshot_dir = Path(*parts[:-1], parts[-1] + '_snaps')
+
 
 class SimulationConfig:
     """Configuration parameters for a standard simulation.
@@ -188,7 +191,7 @@ class SimulationConfig:
                         'cuda': False,
                         'opencl': False,
                         'subgrid': False,
-                        'precision': 'double'}
+                        'precision': 'single'}
 
         self.em_consts = {'c': c, # Speed of light in free space (m/s)
                           'e0': e0, # Permittivity of free space (F/m)
@@ -203,12 +206,14 @@ class SimulationConfig:
             self.general['cuda'] = True
             self.general['cpu'] = False
             self.general['opencl'] = False
+            # Both single and double precision are possible on GPUs, but single
+            # provides best performance.
             self.general['precision'] = 'single'
             self.cuda = {'gpus': [], # gpus: list of GPU objects
                          'gpus_str': [], # gpus_str: list of strings describing GPU(s)
                          'nvcc_opts': None} # nvcc_opts: nvcc compiler options
             # Suppress nvcc warnings on Microsoft Windows
-            if sys.platform == 'win32': self.cuda['nvcc_opts'] = '-w'
+            if sys.platform == 'win32': self.cuda['nvcc_opts'] = ['-w']
 
             # Flatten a list of lists
             if any(isinstance(element, list) for element in self.args.gpu):
@@ -302,7 +307,7 @@ class SimulationConfig:
         # API
         if self.args.inputfile is None:
             self.input_file_path = Path(self.args.outputfile)
-        # CLI
+        # API/CLI
         else:
             self.input_file_path = Path(self.args.inputfile)
 
diff --git a/gprMax/contexts.py b/gprMax/contexts.py
index dfe08bd8..4da8750e 100644
--- a/gprMax/contexts.py
+++ b/gprMax/contexts.py
@@ -72,13 +72,8 @@ class Context:
 
     def print_time_report(self):
         """Print the total simulation time based on context."""
-        s = self.make_time_report()
-        log.info(s)
-
-    def make_time_report(self):
-        """Generate a string for the total simulation time."""
         s = f"\n=== Simulation on {config.sim_config.hostinfo['hostname']} completed in [HH:MM:SS]: {datetime.timedelta(seconds=self.tsimend - self.tsimstart)}"
-        return f"{s} {'=' * (get_terminal_width() - 1 - len(s))}\n"
+        log.info(f"{s} {'=' * (get_terminal_width() - 1 - len(s))}\n")
 
 
 class NoMPIContext(Context):
diff --git a/gprMax/cuda/fields_updates.py b/gprMax/cuda/fields_updates.py
index 6000c016..72cda143 100644
--- a/gprMax/cuda/fields_updates.py
+++ b/gprMax/cuda/fields_updates.py
@@ -160,7 +160,7 @@ __global__ void update_electric_dispersive_A(int NX, int NY, int NZ, int MAXPOLE
         int materialEx = ID[INDEX4D_ID(0,i_ID,j_ID,k_ID)];
         $REAL phi = 0;
         for (int pole = 0; pole < MAXPOLES; pole++) {
-            phi = phi + updatecoeffsdispersive[INDEX2D_MATDISP(materialEx,pole*3)].real() * Tx[INDEX4D_T(pole,i_T,j_T,k_T)].real();
+            phi = phi + updatecoeffsdispersive[INDEX2D_MATDISP(materialEx,pole*3)]$REALFUNC * Tx[INDEX4D_T(pole,i_T,j_T,k_T)]$REALFUNC;
             Tx[INDEX4D_T(pole,i_T,j_T,k_T)] = updatecoeffsdispersive[INDEX2D_MATDISP(materialEx,1+(pole*3))] * Tx[INDEX4D_T(pole,i_T,j_T,k_T)] + updatecoeffsdispersive[INDEX2D_MATDISP(materialEx,2+(pole*3))] * Ex[INDEX3D_FIELDS(i,j,k)];
         }
         Ex[INDEX3D_FIELDS(i,j,k)] = updatecoeffsE[INDEX2D_MAT(materialEx,0)] * Ex[INDEX3D_FIELDS(i,j,k)] + updatecoeffsE[INDEX2D_MAT(materialEx,2)] * (Hz[INDEX3D_FIELDS(i,j,k)] - Hz[INDEX3D_FIELDS(i,j-1,k)]) - updatecoeffsE[INDEX2D_MAT(materialEx,3)] * (Hy[INDEX3D_FIELDS(i,j,k)] - Hy[INDEX3D_FIELDS(i,j,k-1)]) - updatecoeffsE[INDEX2D_MAT(materialEx,4)] * phi;
@@ -171,7 +171,7 @@ __global__ void update_electric_dispersive_A(int NX, int NY, int NZ, int MAXPOLE
         int materialEy = ID[INDEX4D_ID(1,i_ID,j_ID,k_ID)];
         $REAL phi = 0;
         for (int pole = 0; pole < MAXPOLES; pole++) {
-            phi = phi + updatecoeffsdispersive[INDEX2D_MATDISP(materialEy,pole*3)].real() * Ty[INDEX4D_T(pole,i_T,j_T,k_T)].real();
+            phi = phi + updatecoeffsdispersive[INDEX2D_MATDISP(materialEy,pole*3)]$REALFUNC * Ty[INDEX4D_T(pole,i_T,j_T,k_T)]$REALFUNC;
             Ty[INDEX4D_T(pole,i_T,j_T,k_T)] = updatecoeffsdispersive[INDEX2D_MATDISP(materialEy,1+(pole*3))] * Ty[INDEX4D_T(pole,i_T,j_T,k_T)] + updatecoeffsdispersive[INDEX2D_MATDISP(materialEy,2+(pole*3))] * Ey[INDEX3D_FIELDS(i,j,k)];
         }
         Ey[INDEX3D_FIELDS(i,j,k)] = updatecoeffsE[INDEX2D_MAT(materialEy,0)] * Ey[INDEX3D_FIELDS(i,j,k)] + updatecoeffsE[INDEX2D_MAT(materialEy,3)] * (Hx[INDEX3D_FIELDS(i,j,k)] - Hx[INDEX3D_FIELDS(i,j,k-1)]) - updatecoeffsE[INDEX2D_MAT(materialEy,1)] * (Hz[INDEX3D_FIELDS(i,j,k)] - Hz[INDEX3D_FIELDS(i-1,j,k)]) - updatecoeffsE[INDEX2D_MAT(materialEy,4)] * phi;
@@ -182,7 +182,7 @@ __global__ void update_electric_dispersive_A(int NX, int NY, int NZ, int MAXPOLE
         int materialEz = ID[INDEX4D_ID(2,i_ID,j_ID,k_ID)];
         $REAL phi = 0;
         for (int pole = 0; pole < MAXPOLES; pole++) {
-            phi = phi + updatecoeffsdispersive[INDEX2D_MATDISP(materialEz,pole*3)].real() * Tz[INDEX4D_T(pole,i_T,j_T,k_T)].real();
+            phi = phi + updatecoeffsdispersive[INDEX2D_MATDISP(materialEz,pole*3)]$REALFUNC * Tz[INDEX4D_T(pole,i_T,j_T,k_T)]$REALFUNC;
             Tz[INDEX4D_T(pole,i_T,j_T,k_T)] = updatecoeffsdispersive[INDEX2D_MATDISP(materialEz,1+(pole*3))] * Tz[INDEX4D_T(pole,i_T,j_T,k_T)] + updatecoeffsdispersive[INDEX2D_MATDISP(materialEz,2+(pole*3))] * Ez[INDEX3D_FIELDS(i,j,k)];
         }
         Ez[INDEX3D_FIELDS(i,j,k)] = updatecoeffsE[INDEX2D_MAT(materialEz,0)] * Ez[INDEX3D_FIELDS(i,j,k)] + updatecoeffsE[INDEX2D_MAT(materialEz,1)] * (Hy[INDEX3D_FIELDS(i,j,k)] - Hy[INDEX3D_FIELDS(i-1,j,k)]) - updatecoeffsE[INDEX2D_MAT(materialEz,2)] * (Hx[INDEX3D_FIELDS(i,j,k)] - Hx[INDEX3D_FIELDS(i,j-1,k)]) - updatecoeffsE[INDEX2D_MAT(materialEz,4)] * phi;
diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py
index 9a6fe6e0..a26b8d44 100644
--- a/gprMax/model_build_run.py
+++ b/gprMax/model_build_run.py
@@ -76,7 +76,7 @@ class ModelBuildRun:
         self.p = psutil.Process()
 
         # Normal model reading/building process; bypassed if geometry information to be reused
-        self.reuse_geometry() if config.get_model_config().reuse_geometry else self.build_geometry()
+        self.build_geometry() if not config.get_model_config().reuse_geometry else self.reuse_geometry()
 
         log.info(f'\nOutput directory: {config.get_model_config().output_file_path.parent.resolve()}')
 
@@ -136,10 +136,21 @@ class ModelBuildRun:
 
         scene = self.build_scene()
 
-        # Combine available grids and check memory requirements
+        # Combine available grids
         grids = [G] + G.subgrids
-        total_mem = mem_check_all(grids)
-        log.info(f'\nMemory (RAM) required: ~{human_size(total_mem)}')
+
+        # Check for dispersive materials and specific type
+        for grid in grids:
+            if config.get_model_config().materials['maxpoles'] != 0:
+                config.get_model_config().materials['drudelorentz'] = any([m for m in grid.materials if 'drude' in m.type or 'lorentz' in m.type])
+
+        # Set data type if any dispersive materials (must be done before memory checks)
+        if config.get_model_config().materials['maxpoles'] != 0:
+            config.get_model_config().set_dispersive_material_types()
+
+        # Check memory requirements
+        total_mem, mem_strs = mem_check_all(grids)
+        log.info(f'\nMemory required: {" + ".join(mem_strs)} + ~{human_size(config.get_model_config().mem_overhead)} overhead = {human_size(total_mem)}')
 
         # Build grids
         gridbuilders = [GridBuilder(grid) for grid in grids]
@@ -172,7 +183,7 @@ class ModelBuildRun:
     def reuse_geometry(self):
         # Reset iteration number
         self.G.iteration = 0
-        config.get_model_config().set_inputfilestr(f'\n--- Model {config.get_model_config().appendmodelnumber}/{config.sim_config.model_end}, input file (not re-processed, i.e. geometry fixed): {config.sim_config.input_file_path}')
+        config.get_model_config().inputfilestr = f'\n--- Model {config.get_model_config().appendmodelnumber}/{config.sim_config.model_end}, input file (not re-processed, i.e. geometry fixed): {config.sim_config.input_file_path}'
         log.info(config.get_model_config().inputfilestr)
         for grid in [self.G] + self.G.subgrids:
             grid.reset_fields()
@@ -225,11 +236,11 @@ class ModelBuildRun:
             memsolve (float): Memory (RAM) used on GPU.
         """
 
-        memGPU_str = ''
+        mem_str = ''
         if config.sim_config.general['cuda']:
-            memGPU_str = f' host + ~{human_size(memsolve)} GPU'
+            mem_str = f' host + ~{human_size(memsolve)} GPU'
 
-        log.info(f'\nMemory (RAM) used: ~{human_size(self.p.memory_full_info().uss)}{memGPU_str}')
+        log.info(f'\nMemory used: ~{human_size(self.p.memory_full_info().uss)}{mem_str}')
         log.info(f'Solving time [HH:MM:SS]: {datetime.timedelta(seconds=tsolve)}')
 
     def solve(self, solver):
diff --git a/gprMax/solvers.py b/gprMax/solvers.py
index 1de5a38b..690b68d0 100644
--- a/gprMax/solvers.py
+++ b/gprMax/solvers.py
@@ -66,7 +66,6 @@ def create_solver(G):
         updates = CUDAUpdates(G)
         solver = Solver(updates)
 
-
     return solver
 
 
diff --git a/gprMax/subgrids/base.py b/gprMax/subgrids/base.py
index 172af766..b61d88cf 100644
--- a/gprMax/subgrids/base.py
+++ b/gprMax/subgrids/base.py
@@ -25,7 +25,7 @@ from ..grid import FDTDGrid
 class SubGridBase(FDTDGrid):
 
     def __init__(self, *args, **kwargs):
-        super().__init__(args[0])
+        super().__init__()
 
         self.mode = '3D'
         self.ratio = kwargs['ratio']
diff --git a/gprMax/subgrids/subgrid_hsg.py b/gprMax/subgrids/subgrid_hsg.py
index 56a4239c..3befdabd 100644
--- a/gprMax/subgrids/subgrid_hsg.py
+++ b/gprMax/subgrids/subgrid_hsg.py
@@ -16,7 +16,9 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
-from colorama import init, Fore, Style
+from colorama import init
+from colorama import Fore
+from colorama import Style
 init()
 
 import gprMax.config as config
@@ -24,40 +26,43 @@ from .base import SubGridBase
 from ..cython.fields_updates_hsg import cython_update_is
 from ..cython.fields_updates_hsg import cython_update_magnetic_os
 from ..cython.fields_updates_hsg import cython_update_electric_os
-from ..utilities import human_size
 
 
 class SubGridHSG(SubGridBase):
 
     gridtype = '3DSUBGRID'
 
-    def __init__(self, model_num, **kwargs):
-        super().__init__(model_num, **kwargs)
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
         self.gridtype = SubGridHSG.gridtype
 
     def update_magnetic_is(self, precursors):
         """Update the subgrid nodes at the IS with the currents derived
             from the main grid.
-            Args:
-                nwl, nwm, nwn, face, field, inc_field, lookup_id, sign, mod, co
+
+        Args:
+            nwl, nwm, nwn, face, field, inc_field, lookup_id, sign, mod, co
         """
 
         # Hz = c0Hz - c1Ey + c2Ex
         # Hy = c0Hy - c3Ex + c1Ez
         # Hx = c0Hx - c2Ez + c3Ey
         # bottom and top
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx, self.nwy + 1, self.nwz, 1, self.Hy, precursors.ex_bottom, precursors.ex_top, self.IDlookup['Hy'], 1, -1, 3, config.get_model().ompthreads)
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwy, self.nwz, 1, self.Hx, precursors.ey_bottom, precursors.ey_top, self.IDlookup['Hx'], -1, 1, 3, config.get_model().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx, self.nwy + 1, self.nwz, 1, self.Hy, precursors.ex_bottom, precursors.ex_top, self.IDlookup['Hy'], 1, -1, 3, config.get_model_config().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwy, self.nwz, 1, self.Hx, precursors.ey_bottom, precursors.ey_top, self.IDlookup['Hx'], -1, 1, 3, config.get_model_config().ompthreads)
 
         # left and right
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwy, self.nwz + 1, self.nwx, 2, self.Hz, precursors.ey_left, precursors.ey_right, self.IDlookup['Hz'], 1, -1, 1, config.get_model().ompthreads)
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwy + 1, self.nwz, self.nwx, 2, self.Hy, precursors.ez_left, precursors.ez_right, self.IDlookup['Hy'], -1, 1, 1, config.get_model().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwy, self.nwz + 1, self.nwx, 2, self.Hz, precursors.ey_left, precursors.ey_right, self.IDlookup['Hz'], 1, -1, 1, config.get_model_config().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwy + 1, self.nwz, self.nwx, 2, self.Hy, precursors.ez_left, precursors.ez_right, self.IDlookup['Hy'], -1, 1, 1, config.get_model_config().ompthreads)
 
         # front and back
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx, self.nwz + 1, self.nwy, 3, self.Hz, precursors.ex_front, precursors.ex_back, self.IDlookup['Hz'], -1, 1, 2, config.get_model().ompthreads)
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwz, self.nwy, 3, self.Hx, precursors.ez_front, precursors.ez_back, self.IDlookup['Hx'], 1, -1, 2, config.get_model().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx, self.nwz + 1, self.nwy, 3, self.Hz, precursors.ex_front, precursors.ex_back, self.IDlookup['Hz'], -1, 1, 2, config.get_model_config().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsH, self.ID, self.n_boundary_cells, -1, self.nwx + 1, self.nwz, self.nwy, 3, self.Hx, precursors.ez_front, precursors.ez_back, self.IDlookup['Hx'], 1, -1, 2, config.get_model_config().ompthreads)
 
     def update_electric_is(self, precursors):
+        """
+
+        """
         # Args: nwl, nwm, nwn, face, field, inc_field, lookup_id, sign, mod, co
 
         # Ex = c0(Ex) + c2(dHz) - c3(dHy)
@@ -65,18 +70,21 @@ class SubGridHSG(SubGridBase):
         # Ez = c0(Ez) + c1(dHy) - c2(dHx)
 
         # bottom and top
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx, self.nwy + 1, self.nwz, 1, self.Ex, precursors.hy_bottom, precursors.hy_top, self.IDlookup['Ex'], 1, -1, 3, config.get_model().ompthreads)
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx + 1, self.nwy, self.nwz, 1, self.Ey, precursors.hx_bottom, precursors.hx_top, self.IDlookup['Ey'], -1, 1, 3, config.get_model().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx, self.nwy + 1, self.nwz, 1, self.Ex, precursors.hy_bottom, precursors.hy_top, self.IDlookup['Ex'], 1, -1, 3, config.get_model_config().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx + 1, self.nwy, self.nwz, 1, self.Ey, precursors.hx_bottom, precursors.hx_top, self.IDlookup['Ey'], -1, 1, 3, config.get_model_config().ompthreads)
 
         # left and right
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwy, self.nwz + 1, self.nwx, 2, self.Ey, precursors.hz_left, precursors.hz_right, self.IDlookup['Ey'], 1, -1, 1, config.get_model().ompthreads)
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwy + 1, self.nwz, self.nwx, 2, self.Ez, precursors.hy_left, precursors.hy_right, self.IDlookup['Ez'], -1, 1, 1, config.get_model().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwy, self.nwz + 1, self.nwx, 2, self.Ey, precursors.hz_left, precursors.hz_right, self.IDlookup['Ey'], 1, -1, 1, config.get_model_config().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwy + 1, self.nwz, self.nwx, 2, self.Ez, precursors.hy_left, precursors.hy_right, self.IDlookup['Ez'], -1, 1, 1, config.get_model_config().ompthreads)
 
         # front and back
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx, self.nwz + 1, self.nwy, 3, self.Ex, precursors.hz_front, precursors.hz_back, self.IDlookup['Ex'], -1, 1, 2, config.get_model().ompthreads)
-        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx + 1, self.nwz, self.nwy, 3, self.Ez, precursors.hx_front, precursors.hx_back, self.IDlookup['Ez'], 1, -1, 2, config.get_model().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx, self.nwz + 1, self.nwy, 3, self.Ex, precursors.hz_front, precursors.hz_back, self.IDlookup['Ex'], -1, 1, 2, config.get_model_config().ompthreads)
+        cython_update_is(self.nwx, self.nwy, self.nwz, self.updatecoeffsE, self.ID, self.n_boundary_cells, 0, self.nwx + 1, self.nwz, self.nwy, 3, self.Ez, precursors.hx_front, precursors.hx_back, self.IDlookup['Ez'], 1, -1, 2, config.get_model_config().ompthreads)
 
     def update_electric_os(self, main_grid):
+        """
+
+        """
         i_l = self.i0 - self.is_os_sep
         i_u = self.i1 + self.is_os_sep
         j_l = self.j0 - self.is_os_sep
@@ -91,19 +99,21 @@ class SubGridHSG(SubGridBase):
         # Ey = c0(Ey) + c3(dHx) - c1(dHz)
         # Ez = c0(Ez) + c1(dHy) - c2(dHx)
 
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, i_l, i_u, k_l, k_u + 1, j_l, j_u, self.nwy, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hz, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, i_l, i_u + 1, k_l, k_u, j_l, j_u, self.nwy, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hx, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, i_l, i_u, k_l, k_u + 1, j_l, j_u, self.nwy, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hz, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 3, i_l, i_u + 1, k_l, k_u, j_l, j_u, self.nwy, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hx, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
 
         # Left and Right
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, j_l, j_u, k_l, k_u + 1, i_l, i_u, self.nwx, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hz, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, j_l, j_u + 1, k_l, k_u, i_l, i_u, self.nwx, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hy, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, j_l, j_u, k_l, k_u + 1, i_l, i_u, self.nwx, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hz, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 2, j_l, j_u + 1, k_l, k_u, i_l, i_u, self.nwx, main_grid.IDlookup['Ez'], main_grid.Ez, self.Hy, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
 
         # Bottom and Top
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, i_l, i_u, j_l, j_u + 1, k_l, k_u, self.nwz, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hy, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
-        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l, k_u, self.nwz, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hx, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, i_l, i_u, j_l, j_u + 1, k_l, k_u, self.nwz, main_grid.IDlookup['Ex'], main_grid.Ex, self.Hy, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
+        cython_update_electric_os(main_grid.updatecoeffsE, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l, k_u, self.nwz, main_grid.IDlookup['Ey'], main_grid.Ey, self.Hx, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
 
     def update_magnetic_os(self, main_grid):
+        """
 
+        """
         i_l = self.i0 - self.is_os_sep
         i_u = self.i1 + self.is_os_sep
         j_l = self.j0 - self.is_os_sep
@@ -119,28 +129,25 @@ class SubGridHSG(SubGridBase):
         # Args: sub_grid, normal, l_l, l_u, m_l, m_u, n_l, n_u, nwn, lookup_id, field, inc_field, co, sign_n, sign_f):
 
         # Front and back
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, i_l, i_u, k_l, k_u + 1, j_l - 1, j_u, self.nwy, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ex, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, i_l, i_u + 1, k_l, k_u, j_l - 1, j_u, self.nwy, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ez, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, i_l, i_u, k_l, k_u + 1, j_l - 1, j_u, self.nwy, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ex, 2, 1, -1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 3, i_l, i_u + 1, k_l, k_u, j_l - 1, j_u, self.nwy, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ez, 2, -1, 1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
 
         # Left and Right
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, j_l, j_u, k_l, k_u + 1, i_l - 1, i_u, self.nwx, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ey, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, j_l, j_u + 1, k_l, k_u, i_l - 1, i_u, self.nwx, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ez, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, j_l, j_u, k_l, k_u + 1, i_l - 1, i_u, self.nwx, main_grid.IDlookup['Hz'], main_grid.Hz, self.Ey, 1, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 2, j_l, j_u + 1, k_l, k_u, i_l - 1, i_u, self.nwx, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ez, 1, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
 
         # bottom and top
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u, j_l, j_u + 1, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ex, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
-        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ey, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model().ompthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u, j_l, j_u + 1, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hy'], main_grid.Hy, self.Ex, 3, -1, 1, 1, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
+        cython_update_magnetic_os(main_grid.updatecoeffsH, main_grid.ID, 1, i_l, i_u + 1, j_l, j_u, k_l - 1, k_u, self.nwz, main_grid.IDlookup['Hx'], main_grid.Hx, self.Ey, 3, 1, -1, 0, self.ratio, self.is_os_sep, self.n_boundary_cells, config.get_model_config().ompthreads)
 
     def __str__(self):
 
-        mem_use = self.mem_est_basic()
-
         s = '\n'
         s += Fore.CYAN
         s += 'Sub Grid HSG\n'
         s += f'Name: {self.name}\n'
-        s += f'dx, dy, dz: {self.dx}m {self.dy}m {self.dz}m\n'
-        s += f'dt: {self.dt}s\n'
-        s += f'Memory Estimate: {human_size(mem_use)}\n'
+        s += f'Spatial discretisation: {self.dx:g} x {self.dy:g} x {self.dz:g}m\n'
+        s += f'Time step (at CFL limit): {self.dt:g} secs\n'
         s += f'Position: ({self.x1}m, {self.y1}m, {self.z1}m), ({self.x2}m, {self.y2}m, {self.z2}m)\n'
         s += f'Main Grid Indices: lower left({self.i0}, {self.j0}, {self.k0}), upper right({self.i1}, {self.j1}, {self.k1})\n'
         s += f'Total Cells: {self.nx} {self.ny} {self.nz}\n'
diff --git a/gprMax/subgrids/updates.py b/gprMax/subgrids/updates.py
index 5185d07d..24147955 100644
--- a/gprMax/subgrids/updates.py
+++ b/gprMax/subgrids/updates.py
@@ -16,19 +16,23 @@
 # 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
+
 from ..exceptions import GeneralError
 from .precursor_nodes import PrecursorNodes
 from .precursor_nodes import PrecursorNodesFiltered
 from .subgrid_hsg import SubGridHSG
 from ..updates import CPUUpdates
 
+log = logging.getLogger(__name__)
+
 
 def create_updates(G):
     """Return the solver for the given subgrids."""
     updaters = []
 
     for sg in G.subgrids:
-        print(sg)
+        log.info(sg)
         sg_type = type(sg)
         if sg_type == SubGridHSG and sg.filter:
             precursors = PrecursorNodesFiltered(G, sg)
diff --git a/gprMax/subgrids/user_objects.py b/gprMax/subgrids/user_objects.py
index 7421f959..fc7eaba6 100644
--- a/gprMax/subgrids/user_objects.py
+++ b/gprMax/subgrids/user_objects.py
@@ -169,7 +169,6 @@ class SubGridHSG(SubGridBase):
                  interpolation=1,
                  filter=True,
                  **kwargs):
-        """Constructor."""
 
         pml_separation = ratio // 2 + 2
 
@@ -189,7 +188,7 @@ class SubGridHSG(SubGridBase):
         self.hash = '#subgrid_hsg'
 
     def create(self, grid, uip):
-        sg = SubGridHSGUser(grid.model_num, **self.kwargs)
+        sg = SubGridHSGUser(**self.kwargs)
         self.setup(sg, grid, uip)
         return sg
 
diff --git a/gprMax/updates.py b/gprMax/updates.py
index 182ce8ce..93a5564f 100644
--- a/gprMax/updates.py
+++ b/gprMax/updates.py
@@ -279,7 +279,8 @@ class CUDAUpdates:
         if config.get_model_config().materials['maxpoles'] > 0:
             kernels_fields = self.source_module(kernel_template_fields.substitute(
                                                 REAL=config.sim_config.dtypes['C_float_or_double'],
-                                                COMPLEX=config.sim_config.dtypes['C_complex'],
+                                                REALFUNC=config.get_model_config().materials['cuda_real_func'],
+                                                COMPLEX=config.get_model_config().materials['dispersiveCdtype'],
                                                 N_updatecoeffsE=self.grid.updatecoeffsE.size,
                                                 N_updatecoeffsH=self.grid.updatecoeffsH.size,
                                                 NY_MATCOEFFS=self.grid.updatecoeffsE.shape[1],
@@ -294,9 +295,11 @@ class CUDAUpdates:
                                                 NY_T=self.grid.Tx.shape[2],
                                                 NZ_T=self.grid.Tx.shape[3]),
                                                 options=config.sim_config.cuda['nvcc_opts'])
-        else: # Set to one any substitutions for dispersive materials
+        else: # Set to one any substitutions for dispersive materials.
+              # Value of COMPLEX is not relevant.
             kernels_fields = self.source_module(kernel_template_fields.substitute(
                                                 REAL=config.sim_config.dtypes['C_float_or_double'],
+                                                REALFUNC=config.get_model_config().materials['cuda_real_func'],
                                                 COMPLEX=config.sim_config.dtypes['C_complex'],
                                                 N_updatecoeffsE=self.grid.updatecoeffsE.size,
                                                 N_updatecoeffsH=self.grid.updatecoeffsH.size,
diff --git a/gprMax/user_inputs.py b/gprMax/user_inputs.py
index 578da69c..2d4e6b88 100644
--- a/gprMax/user_inputs.py
+++ b/gprMax/user_inputs.py
@@ -73,11 +73,9 @@ class UserInput:
             # Incorrect index
             i = p[v.index(err.args[0])]
             if name:
-                s = f"'{cmd_str}' the {err.args[0]} {name}-coordinate {i * dl:g} \
-                    is not within the model domain"
+                s = f"'{cmd_str}' the {err.args[0]} {name}-coordinate {i * dl:g} is not within the model domain"
             else:
-                s = f"'{cmd_str}' {err.args[0]}-coordinate {i * dl:g} is not \
-                    within the model domain"
+                s = f"'{cmd_str}' {err.args[0]}-coordinate {i * dl:g} is not within the model domain"
             raise CmdInputError(log.exception(s))
 
     def discretise_point(self, p):
@@ -106,8 +104,7 @@ class MainGridUserInput(UserInput):
         p = self.check_point(p, cmd_str, name)
 
         if self.grid.within_pml(p):
-            log.warning(Fore.RED + f"'{cmd_str}' sources and receivers should \
-                        not normally be positioned within the PML." + Style.RESET_ALL)
+            log.warning(Fore.RED + f"'{cmd_str}' sources and receivers should not normally be positioned within the PML." + Style.RESET_ALL)
 
         return p
 
@@ -116,8 +113,7 @@ class MainGridUserInput(UserInput):
         p2 = self.check_point(p2, cmd_str, name='upper')
 
         if np.greater(p1, p2).any():
-            raise CmdInputError(log.exception(f"'{cmd_str}' the lower coordinates \
-                                should be less than the upper coordinates."))
+            raise CmdInputError(log.exception(f"'{cmd_str}' the lower coordinates should be less than the upper coordinates."))
 
         return p1, p2
 
@@ -139,10 +135,9 @@ class SubgridUserInput(MainGridUserInput):
         super().__init__(grid)
 
         # Defines the region exposed to the user
-        self.inner_bound = np.array([
-                                    grid.n_boundary_cells_x,
-                                    grid.n_boundary_cells_y,
-                                    grid.n_boundary_cells_z])
+        self.inner_bound = np.array([grid.n_boundary_cells_x,
+                                     grid.n_boundary_cells_y,
+                                     grid.n_boundary_cells_z])
 
         self.outer_bound = np.subtract([grid.nx, grid.ny, grid.nz],
                                        self.inner_bound)
@@ -180,6 +175,5 @@ class SubgridUserInput(MainGridUserInput):
         # the OS non-working region.
         if (np.less(p_t, self.inner_bound).any() or
             np.greater(p_t, self.outer_bound).any()):
-                log.warning(Fore.RED + f"'{cmd_str}' this object traverses the \
-                            Outer Surface. This is an advanced feature." + Style.RESET_ALL)
+                log.warning(Fore.RED + f"'{cmd_str}' this object traverses the Outer Surface. This is an advanced feature." + Style.RESET_ALL)
         return p_t
diff --git a/gprMax/utilities.py b/gprMax/utilities.py
index e68892cb..12758aac 100644
--- a/gprMax/utilities.py
+++ b/gprMax/utilities.py
@@ -28,7 +28,12 @@ import subprocess
 from shutil import get_terminal_size
 import sys
 import textwrap
-from time import perf_counter
+
+try:
+    from time import thread_time as timer_fn
+except ImportError:
+    from time import perf_counter as timer_fn
+
 
 from colorama import init
 from colorama import Fore
@@ -301,7 +306,7 @@ def get_host_info():
             cpuIDinfo = subprocess.check_output("cat /proc/cpuinfo", shell=True, stderr=subprocess.STDOUT, env=myenv).decode('utf-8').strip()
             for line in cpuIDinfo.split('\n'):
                 if re.search('model name', line):
-                    cpuID = re.sub('.*model name.*:', '', line, 1).strip()            
+                    cpuID = re.sub('.*model name.*:', '', line, 1).strip()
             allcpuinfo = subprocess.check_output("lscpu", shell=True, stderr=subprocess.STDOUT, env=myenv).decode('utf-8').strip()
             for line in allcpuinfo.split('\n'):
                 if 'Socket(s)' in line:
@@ -420,44 +425,44 @@ def mem_check_all(grids):
 
     Returns:
         total_mem (int): Total memory required for all grids.
+        mem_strs (list): Strings containing text of memory requirements for
+                            each grid.
     """
 
-    total_mem = 0
     total_snaps_mem = 0
+    mem_strs = []
 
     for grid in grids:
+        # Memory required for main grid arrays
         config.get_model_config().mem_use += grid.mem_est_basic()
+        grid.mem_use += grid.mem_est_basic()
 
-        # Set datatype for dispersive arrays if there are any dispersive materials.
+        # Additional memory required if there are any dispersive materials.
         if config.get_model_config().materials['maxpoles'] != 0:
-            drudelorentz = any([m for m in grid.materials if 'drude' in m.type or 'lorentz' in m.type])
-            if drudelorentz:
-                config.get_model_config().materials['dispersivedtype'] = config.sim_config.dtypes['complex']
-                config.get_model_config().materials['dispersiveCdtype'] = config.sim_config.dtypes['C_complex']
-            else:
-                config.get_model_config().materials['dispersivedtype'] = config.sim_config.dtypes['float_or_double']
-                config.get_model_config().materials['dispersiveCdtype'] = config.sim_config.dtypes['C_float_or_double']
-
-            # Update estimated memory (RAM) usage
             config.get_model_config().mem_use += grid.mem_est_dispersive()
+            grid.mem_use += grid.mem_est_dispersive()
 
-        # Calculate snapshot memory
+        # Additional memory required if there are any snapshots
         if grid.snapshots:
             for snap in grid.snapshots:
                 # 2 x required to account for electric and magnetic fields
-                config.get_model_config().mem_use += int(2 * snap.datasizefield)
-                total_snaps_mem += int(2 * snap.datasizefield)
+                snap_mem = int(2 * snap.datasizefield)
+                config.get_model_config().mem_use += snap_mem
+                total_snaps_mem += snap_mem
+                grid.mem_use += snap_mem
 
-        total_mem += config.get_model_config().mem_use
+        mem_strs.append(f'~{human_size(grid.mem_use)} ({grid.name})')
+
+    total_mem = config.get_model_config().mem_use
 
     # Check if there is sufficient memory on host
     mem_check_host(total_mem)
 
     # Check if there is sufficient memory for any snapshots on GPU
-    if config.sim_config.general['cuda']:
+    if total_snaps_mem > 0 and config.sim_config.general['cuda']:
         mem_check_gpu_snaps(total_mem, total_snaps_mem)
 
-    return total_mem
+    return total_mem, mem_strs
 
 
 class GPU:
@@ -535,5 +540,6 @@ def detect_check_gpus(deviceIDs):
 
 
 def timer():
-    """Function to return the current process wide time in fractional seconds."""
-    return perf_counter()
+    """Function to return time in fractional seconds."""
+    log.debug('Review "thread_time" not currently available in macOS and bug (https://bugs.python.org/issue36205) with "process_time"')
+    return timer_fn()