Basic cylinder_Ascan_2D model is running, but still work to do.

gprMax/cmds_single_use.py

@@ -20,7 +20,6 @@ import decimal as d
 import inspect
 import logging
 import os
-import sys
 
 from colorama import init
 from colorama import Fore
@@ -34,7 +33,6 @@ import gprMax.config as config
 from .exceptions import CmdInputError
 from .waveforms import Waveform
 from .utilities import round_value
-from .utilities import set_openmp_threads
 
 log = logging.getLogger(__name__)
 
@@ -118,15 +116,6 @@ class Domain(UserObjectSingle):
         log.info(f'Mode: {G.mode}')
         log.info(f'Time step (at CFL limit): {G.dt:g} secs')
 
-        # Number of threads (OpenMP) to use
-        G.nthreads = set_openmp_threads()
-        log.info(f'Number of CPU (OpenMP) threads: {G.nthreads}')
-
-        if G.nthreads > config.sim_config.hostinfo['physicalcores']:
-            log.warning(Fore.RED + f"You have specified more threads ({G.nthreads}) \
-                        than available physical CPU cores ({config.hostinfo['physicalcores']}). \
-                        This may lead to degraded performance." + Style.RESET_ALL)
-
 
 class Discretisation(UserObjectSingle):
     """Allows you to specify the discretization of space in the x, y, and z
@@ -288,29 +277,14 @@ class NumThreads(UserObjectSingle):
             return '#n_threads:'
 
     def create(self, G, uip):
-        # Get information about host machine
-
         try:
             n = self.kwargs['n']
         except KeyError:
             raise CmdInputError(self.__str__() + ' requires exactly one parameter to specify the number of threads to use')
-
         if n < 1:
             raise CmdInputError(self.__str__() + ' requires the value to be an integer not less than one')
 
-        G.nthreads = n
-        os.environ['OMP_NUM_THREADS'] = str(G.nthreads)
-
-        log.info(f'Number of CPU (OpenMP) threads: {G.nthreads}')
-        if G.nthreads > config.hostinfo['physicalcores']:
-            log.warning(Fore.RED + f"You have specified more threads ({G.nthreads}) \
-                        than available physical CPU cores (\
-                        {config.sim_config.hostinfo['physicalcores']}). \
-                        This may lead to degraded performance." + Style.RESET_ALL)
-
-        # Print information about any GPU in use
-        if G.gpu is not None:
-            log.info(f'GPU solving using: {G.gpu.deviceID} - {G.gpu.name}')
+        config.model_configs[G.model_num].ompthreads = set_omp_threads(n)
 
 
 class TimeStepStabilityFactor(UserObjectSingle):

gprMax/config.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
 from pathlib import Path
 
 from colorama import init
@@ -31,6 +32,9 @@ from scipy.constants import mu_0 as m0
 from .utilities import get_host_info
 from .utilities import get_terminal_width
 
+log = logging.getLogger(__name__)
+
+
 # Single instance of SimConfig to hold simulation configuration parameters.
 sim_config = None
 
@@ -42,13 +46,22 @@ class ModelConfig:
         N.B. Multiple models can exist within a simulation
     """
 
-    def __init__(self, i):
+    def __init__(self, model_num):
         """
         Args:
-            i (int): Current model number.
+            model_num (int): Model number.
         """
 
-        self.i = i # Indexed from 0
+        self.i = model_num # Indexed from 0
+        self.ompthreads = None # Number of OpenMP threads
+
+        # Store information for CUDA solver type
+        #   gpu: GPU object
+        #   snapsgpu2cpu: copy snapshot data from GPU to CPU during simulation
+        #     N.B. This will happen if the requested snapshots are too large to fit
+        #     on the memory of the GPU. If True this will slow performance significantly
+        self.cuda = {'gpu': None, 'snapsgpu2cpu': False}
+        self.memoryusage = 0 # Total memory usage for all grids in the model
         self.reuse_geometry = False
 
         if not sim_config.single_model:
@@ -65,9 +78,8 @@ class ModelConfig:
         self.snapshot_dir = '_snaps'
 
         # String to print at start of each model run
-        inputfilestr_f = '\n--- Model {}/{}, input file: {}'
-        self.inputfilestr = inputfilestr_f.format(self.i + 1, sim_config.model_end, sim_config.input_file_path)
-        self.next_model = Fore.GREEN + '{} {}\n'.format(self.inputfilestr, '-' * (get_terminal_width() - 1 - len(self.inputfilestr))) + Style.RESET_ALL
+        inputfilestr = f'\n--- Model {self.i + 1}/{sim_config.model_end}, input file: {sim_config.input_file_path}'
+        self.next_model = Fore.GREEN + f"{inputfilestr} {'-' * (get_terminal_width() - 1 - len(inputfilestr))}\n" + Style.RESET_ALL
 
         # Numerical dispersion analysis parameters
         #   highestfreqthres: threshold (dB) down from maximum power (0dB) of main frequency used
@@ -86,6 +98,13 @@ class ModelConfig:
                           'dispersivedtype': None,
                           'dispersiveCdtype': None}
 
+    def memory_check(self):
+        """Check if the required amount of memory (RAM) is available to build
+            and/or run model on the host.
+        """
+        if self.memoryusage > config.sim_config.hostinfo['ram']:
+            raise GeneralError(f"Memory (RAM) required ~{human_size(self.memoryusage)} exceeds {human_size(config.hostinfo['ram'], a_kilobyte_is_1024_bytes=True)} detected!\n")
+
     def get_scene(self):
         if sim_config.scenes:
             return sim_config.scenes[self.i]
@@ -113,6 +132,7 @@ class SimulationConfig:
         """
 
         self.args = args
+        log.debug('Fix parsing args')
 
         # General settings for the simulation
         #   inputfilepath: path to inputfile location
@@ -143,12 +163,8 @@ class SimulationConfig:
         # Store information about host machine
         self.hostinfo = get_host_info()
 
-        # Store information for CUDA solver type
-        #   gpus: information about any GPUs as a list of GPU objects
-        #   snapsgpu2cpu: copy snapshot data from GPU to CPU during simulation
-        #     N.B. This will happen if the requested snapshots are too large to fit
-        #     on the memory of the GPU. If True this will slow performance significantly
-        self.cuda = {'gpus': None, 'snapsgpu2cpu': False}
+        # Information about any GPUs as a list of GPU objects
+        self.cuda_gpus = []
 
         # Data type (precision) for electromagnetic field output
         self.dtypes = None

gprMax/contexts.py

@@ -17,6 +17,7 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 import datetime
+import logging
 
 import gprMax.config as config
 from ._version import __version__, codename
@@ -25,9 +26,12 @@ from .model_build_run import ModelBuildRun
 from .solvers import create_solver
 from .solvers import create_G
 from .utilities import get_terminal_width
+from .utilities import human_size
 from .utilities import logo
 from .utilities import timer
 
+log = logging.getLogger(__name__)
+
 
 class Context:
     """Generic context for the model to run in. Sub-class with specific contexts
@@ -42,14 +46,21 @@ class Context:
     def run(self):
         """Run the simulation in the correct context."""
         self.print_logo_copyright()
+        self.print_hostinfo()
         self.tsimstart = timer()
         self._run()
         self.tsimend = timer()
+        self.print_time_report()
 
     def print_logo_copyright(self):
         """Print gprMax logo, version, and copyright/licencing information."""
         logo(__version__ + ' (' + codename + ')')
 
+    def print_hostinfo(self):
+        """Print information about the host machine."""
+        hyperthreadingstr = f", {config.sim_config.hostinfo['logicalcores']} cores with Hyper-Threading" if config.sim_config.hostinfo['hyperthreading'] else ''
+        log.info(f"\nHost: {config.sim_config.hostinfo['hostname']} | {config.sim_config.hostinfo['machineID']} | {config.sim_config.hostinfo['sockets']} x {config.sim_config.hostinfo['cpuID']} ({config.sim_config.hostinfo['physicalcores']} cores{hyperthreadingstr}) | {human_size(config.sim_config.hostinfo['ram'], a_kilobyte_is_1024_bytes=True)} RAM | {config.sim_config.hostinfo['osversion']}")
+
     def print_time_report(self):
         """Print the total simulation time based on context."""
         s = self.make_time_report()
@@ -57,7 +68,8 @@ class Context:
 
     def make_time_report(self):
         """Generate a string for the total simulation time."""
-        pass
+        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"
 
 
 class NoMPIContext(Context):
@@ -75,26 +87,17 @@ class NoMPIContext(Context):
             # The next model to run only gets a new solver if the
             # geometry is not re-used.
             if i != 0 and config.sim_config.args.geometry_fixed:
-                config.model_config[i].reuse_geometry = True
+                config.model_configs[i].reuse_geometry = True
             else:
-                G = create_G()
+                G = create_G(i)
 
-            model = ModelBuildRun(i, G)
+            model = ModelBuildRun(G)
             model.build()
             solver = create_solver(G)
 
             if not config.sim_config.args.geometry_only:
                 model.solve(solver)
 
-    def make_time_report(self):
-        """Function to specialise the time reporting for the standard simulation
-            context.
-        """
-
-        sim_time = datetime.timedelta(seconds=self.tsimend - self.tsimstart)
-        s = f"\n=== Simulation on {self.simconfig.hostinfo['hostname']} completed in [HH:MM:SS]: {sim_time}"
-        return f"{s} {'=' * (get_terminal_width() - 1 - len(s))}\n"
-
 
 class MPIContext(Context):
     """Mixed mode MPI/OpenMP/CUDA context - MPI task farm is used to distribute
@@ -109,9 +112,6 @@ class MPIContext(Context):
     def _run(self):
         pass
 
-    def make_time_report(self):
-        pass
-
 
 class MPINoSpawnContext(Context):
 
@@ -122,9 +122,6 @@ class MPINoSpawnContext(Context):
     def _run(self):
         pass
 
-    def make_time_report(self):
-        pass
-
 
 def create_context():
     """Create a context in which to run the simulation. i.e MPI.

gprMax/geometry_outputs.py

@@ -17,6 +17,7 @@
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
 import os
+from pathlib import Path
 import sys
 
 import h5py
@@ -114,8 +115,9 @@ class GeometryView:
             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
+        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):
         """

gprMax/gprMax.py

@@ -140,4 +140,3 @@ def run_main(args):
     write_simulation_config(args)
     context = create_context()
     context.run()
-    context.print_time_report()

gprMax/grid.py

@@ -41,13 +41,15 @@ class FDTDGrid:
         accessing regularly used parameters.
     """
 
-    def __init__(self):
+    def __init__(self, model_num):
+        """
+        Args:
+            model_num (int): Model number.
+        """
+
         self.title = ''
-        self.memoryusage = 0
         self.name = 'Main'
-        self.mode = '' # 2D TMx, 2D TMy, 2D TMz, or 3D
-        self.gpu = None
-        self.outputdirectory = ''
+        self.model_num = model_num
 
         self.nx = 0
         self.ny = 0
@@ -206,14 +208,6 @@ class FDTDGrid:
 
         self.memoryusage = int(stdoverhead + fieldarrays + solidarray + rigidarrays + pmlarrays)
 
-    def memory_check(self):
-        """Check if the required amount of memory (RAM) is available to build
-            and/or run model on the host.
-        """
-        if self.memoryusage > config.sim_config.hostinfo['ram']:
-            raise GeneralError(f"Memory (RAM) required ~{human_size(self.memoryusage)} \
-                               exceeds {human_size(config.hostinfo['ram'], a_kilobyte_is_1024_bytes=True)} detected!\n")
-
     def tmx(self):
         """Add PEC boundaries to invariant direction in 2D TMx mode.
             N.B. 2D modes are a single cell slice of 3D grid.
@@ -370,7 +364,7 @@ def dispersion_analysis(G):
 
                     # Set maximum frequency to a threshold drop from maximum power, ignoring DC value
                     try:
-                        freqthres = np.where(power[freqmaxpower:] < -config.numdispersion['highestfreqthres'])[0][0] + freqmaxpower
+                        freqthres = np.where(power[freqmaxpower:] < -config.model_configs[G.model_num].numdispersion['highestfreqthres'])[0][0] + freqmaxpower
                         results['maxfreq'].append(freqs[freqthres])
                     except ValueError:
                         results['error'] = 'unable to calculate maximum power from waveform, most likely due to undersampling.'
@@ -411,9 +405,9 @@ def dispersion_analysis(G):
         minwavelength = minvelocity / results['maxfreq']
 
         # Maximum spatial step
-        if '3D' in config.general['mode']:
+        if '3D' in config.sim_config.general['mode']:
             delta = max(G.dx, G.dy, G.dz)
-        elif '2D' in config.general['mode']:
+        elif '2D' in config.sim_config.general['mode']:
             if G.nx == 1:
                 delta = max(G.dy, G.dz)
             elif G.ny == 1:
@@ -428,7 +422,7 @@ def dispersion_analysis(G):
         results['N'] = minwavelength / delta
 
         # Check grid sampling will result in physical wave propagation
-        if int(np.floor(results['N'])) >= config.numdispersion['mingridsampling']:
+        if int(np.floor(results['N'])) >= config.model_configs[G.model_num].numdispersion['mingridsampling']:
             # Numerical phase velocity
             vp = np.pi / (results['N'] * np.arcsin((1 / S) * np.sin((np.pi * S) / results['N'])))
 

gprMax/materials.py

@@ -131,12 +131,18 @@ class DispersiveMaterial(Material):
 
         # The implementation of the dispersive material modelling comes from the
         # derivation in: http://dx.doi.org/10.1109/TAP.2014.2308549
-        self.w = np.zeros(config.materials['maxpoles'], dtype=config.materials['dispersivedtype'])
-        self.q = np.zeros(config.materials['maxpoles'], dtype=config.materials['dispersivedtype'])
-        self.zt = np.zeros(config.materials['maxpoles'], dtype=config.materials['dispersivedtype'])
-        self.zt2 = np.zeros(config.materials['maxpoles'], dtype=config.materials['dispersivedtype'])
-        self.eqt = np.zeros(config.materials['maxpoles'], dtype=config.materials['dispersivedtype'])
-        self.eqt2 = np.zeros(config.materials['maxpoles'], dtype=config.materials['dispersivedtype'])
+        self.w = np.zeros(config.model_configs[G.model_num].materials['maxpoles'],
+                          dtype=config.materials['dispersivedtype'])
+        self.q = np.zeros(config.model_configs[G.model_num].materials['maxpoles'],
+                          dtype=config.materials['dispersivedtype'])
+        self.zt = np.zeros(config.model_configs[G.model_num].materials['maxpoles'],
+                           dtype=config.materials['dispersivedtype'])
+        self.zt2 = np.zeros(config.model_configs[G.model_num].materials['maxpoles'],
+                            dtype=config.materials['dispersivedtype'])
+        self.eqt = np.zeros(config.model_configs[G.model_num].materials['maxpoles'],
+                            dtype=config.materials['dispersivedtype'])
+        self.eqt2 = np.zeros(config.model_configs[G.model_num].materials['maxpoles'],
+                            dtype=config.materials['dispersivedtype'])
 
         for x in range(self.poles):
             if 'debye' in self.type:
@@ -217,7 +223,7 @@ def process_materials(G):
                                 print a table.
     """
 
-    if config.materials['maxpoles'] == 0:
+    if config.model_configs[G.model_num].materials['maxpoles'] == 0:
         materialsdata = [['\nID', '\nName', '\nType', '\neps_r', 'sigma\n[S/m]',
                           '\nmu_r', 'sigma*\n[Ohm/m]', 'Dielectric\nsmoothable']]
     else:
@@ -237,7 +243,7 @@ def process_materials(G):
         # Add update coefficients to overall storage for dispersive materials
         if hasattr(material, 'poles'):
             z = 0
-            for pole in range(config.materials['maxpoles']):
+            for pole in range(config.model_configs[G.model_num].materials['maxpoles']):
                 G.updatecoeffsdispersive[material.numID, z:z + 3] = (config.sim_config.em_consts['e0'] *
                                                                      material.eqt2[pole], material.eqt[pole], material.zt[pole])
                 z += 3
@@ -249,7 +255,7 @@ def process_materials(G):
         materialtext.append(material.type)
         materialtext.append('{:g}'.format(material.er))
         materialtext.append('{:g}'.format(material.se))
-        if config.materials['maxpoles'] > 0:
+        if config.model_configs[G.model_num].materials['maxpoles'] > 0:
             if 'debye' in material.type:
                 materialtext.append('\n'.join('{:g}'.format(deltaer) for deltaer in material.deltaer))
                 materialtext.append('\n'.join('{:g}'.format(tau) for tau in material.tau))
@@ -361,8 +367,8 @@ class PeplinskiSoil:
                 m.type = 'debye'
                 m.averagable = False
                 m.poles = 1
-                if m.poles > config.materials['maxpoles']:
-                    config.materials['maxpoles'] = m.poles
+                if m.poles > config.model_configs[G.model_num].materials['maxpoles']:
+                    config.model_configs[G.model_num].materials['maxpoles'] = m.poles
                 m.er = eri
                 m.se = sig
                 m.deltaer.append(er - eri)

gprMax/model_build_run.py

@@ -71,6 +71,7 @@ from .utilities import get_terminal_width
 from .utilities import human_size
 from .utilities import open_path_file
 from .utilities import round32
+from .utilities import set_omp_threads
 from .utilities import timer
 
 log = logging.getLogger(__name__)
@@ -79,8 +80,7 @@ log = logging.getLogger(__name__)
 class ModelBuildRun:
     """Builds and runs (solves) a model."""
 
-    def __init__(self, i, G):
-        self.i = i
+    def __init__(self, G):
         self.G = G
         # Monitor memory usage
         self.p = None
@@ -88,55 +88,67 @@ class ModelBuildRun:
     def build(self):
         """Builds the Yee cells for a model."""
 
+        G = self.G
+
         # Monitor memory usage
         self.p = psutil.Process()
 
         # Normal model reading/building process; bypassed if geometry information to be reused
-        self.reuse_geometry() if config.model_configs[self.i].reuse_geometry else self.build_geometry()
-
-        G = self.G
+        self.reuse_geometry() if config.model_configs[G.model_num].reuse_geometry else self.build_geometry()
 
         # Adjust position of simple sources and receivers if required
         if G.srcsteps[0] != 0 or G.srcsteps[1] != 0 or G.srcsteps[2] != 0:
             for source in itertools.chain(G.hertziandipoles, G.magneticdipoles):
-                if config.model_configs[self.i].i == 0:
-                    if source.xcoord + G.srcsteps[0] * self.sim_config.model_end < 0 or source.xcoord + G.srcsteps[0] * self.sim_config.model_end > G.nx or source.ycoord + G.srcsteps[1] * self.sim_config.model_end < 0 or source.ycoord + G.srcsteps[1] * self.sim_config.model_end > G.ny or source.zcoord + G.srcsteps[2] * self.sim_config.model_end < 0 or source.zcoord + G.srcsteps[2] * self.sim_config.model_end > G.nz:
+                if config.model_configs[G.model_num] == 0:
+                    if (source.xcoord + G.srcsteps[0] * self.sim_config.model_end < 0 or
+                        source.xcoord + G.srcsteps[0] * self.sim_config.model_end > G.nx or
+                        source.ycoord + G.srcsteps[1] * self.sim_config.model_end < 0 or
+                        source.ycoord + G.srcsteps[1] * self.sim_config.model_end > G.ny or
+                        source.zcoord + G.srcsteps[2] * self.sim_config.model_end < 0 or
+                        source.zcoord + G.srcsteps[2] * self.sim_config.model_end > G.nz):
                         raise GeneralError('Source(s) will be stepped to a position outside the domain.')
-                source.xcoord = source.xcoordorigin + (config.model_configs[self.i].i) * G.srcsteps[0]
-                source.ycoord = source.ycoordorigin + (config.model_configs[self.i].i) * G.srcsteps[1]
-                source.zcoord = source.zcoordorigin + (config.model_configs[self.i].i) * G.srcsteps[2]
+                source.xcoord = source.xcoordorigin + config.model_configs[G.model_num] * G.srcsteps[0]
+                source.ycoord = source.ycoordorigin + config.model_configs[G.model_num] * G.srcsteps[1]
+                source.zcoord = source.zcoordorigin + config.model_configs[G.model_num] * G.srcsteps[2]
         if G.rxsteps[0] != 0 or G.rxsteps[1] != 0 or G.rxsteps[2] != 0:
             for receiver in G.rxs:
-                if config.model_configs[self.i].i == 0:
-                    if receiver.xcoord + G.rxsteps[0] * self.sim_config.model_end < 0 or receiver.xcoord + G.rxsteps[0] * self.sim_config.model_end > G.nx or receiver.ycoord + G.rxsteps[1] * self.sim_config.model_end < 0 or receiver.ycoord + G.rxsteps[1] * self.sim_config.model_end > G.ny or receiver.zcoord + G.rxsteps[2] * self.sim_config.model_end < 0 or receiver.zcoord + G.rxsteps[2] * self.sim_config.model_end > G.nz:
+                if config.model_configs[G.model_num] == 0:
+                    if (receiver.xcoord + G.rxsteps[0] * self.sim_config.model_end < 0 or
+                        receiver.xcoord + G.rxsteps[0] * self.sim_config.model_end > G.nx or
+                        receiver.ycoord + G.rxsteps[1] * self.sim_config.model_end < 0 or
+                        receiver.ycoord + G.rxsteps[1] * self.sim_config.model_end > G.ny or
+                        receiver.zcoord + G.rxsteps[2] * self.sim_config.model_end < 0 or
+                        receiver.zcoord + G.rxsteps[2] * self.sim_config.model_end > G.nz):
                         raise GeneralError('Receiver(s) will be stepped to a position outside the domain.')
-                receiver.xcoord = receiver.xcoordorigin + (config.model_configs[self.i].i) * G.rxsteps[0]
-                receiver.ycoord = receiver.ycoordorigin + (config.model_configs[self.i].i) * G.rxsteps[1]
-                receiver.zcoord = receiver.zcoordorigin + (config.model_configs[self.i].i) * G.rxsteps[2]
+                receiver.xcoord = receiver.xcoordorigin + config.model_configs[G.model_num] * G.rxsteps[0]
+                receiver.ycoord = receiver.ycoordorigin + config.model_configs[G.model_num] * G.rxsteps[1]
+                receiver.zcoord = receiver.zcoordorigin + config.model_configs[G.model_num] * G.rxsteps[2]
 
         # Write files for any geometry views and geometry object outputs
         if not (G.geometryviews or G.geometryobjectswrite) and self.sim_config.geometry_only and config.is_messages():
             log.warning(Fore.RED + f'\nNo geometry views or geometry objects to output found.' + Style.RESET_ALL)
-        if config.is_messages(): print()
+        if config.sim_config.is_messages(): log.info('')
         for i, geometryview in enumerate(G.geometryviews):
-            geometryview.set_filename(config.model_configs[self.i].appendmodelnumber)
-            pbar = tqdm(total=geometryview.datawritesize, unit='byte', unit_scale=True, desc='Writing geometry view file {}/{}, {}'.format(i + 1, len(G.geometryviews), os.path.split(geometryview.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars'])
+            geometryview.set_filename(config.model_configs[G.model_num].appendmodelnumber)
+            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,
+                        disable=not config.sim_config.general['progressbars'])
             geometryview.write_vtk(G, pbar)
             pbar.close()
         for i, geometryobject in enumerate(G.geometryobjectswrite):
-            pbar = tqdm(total=geometryobject.datawritesize, unit='byte', unit_scale=True, desc='Writing geometry object file {}/{}, {}'.format(i + 1, len(G.geometryobjectswrite), os.path.split(geometryobject.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars'])
+            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,
+                        disable=not config.sim_config.general['progressbars'])
             geometryobject.write_hdf5(G, pbar)
             pbar.close()
-        if config.is_messages(): print()
-
-        # # If only writing geometry information
-        # if self.sim_config.geometry_only:
-        #     tsolve = 0
+        if config.sim_config.is_messages(): log.info('')
 
     def build_geometry(self):
         G = self.G
 
-        log.info(config.model_configs[self.i].next_model)
+        log.info(config.model_configs[G.model_num].next_model)
 
         scene = self.build_scene()
 
@@ -153,7 +165,7 @@ class ModelBuildRun:
             gb.grid.initialise_std_update_coeff_arrays()
 
         # Set datatype for dispersive arrays if there are any dispersive materials.
-        if config.model_configs[self.i].materials['maxpoles'] != 0:
+        if config.model_configs[G.model_num].materials['maxpoles'] != 0:
             drudelorentz = any([m for m in G.materials if 'drude' in m.type or 'lorentz' in m.type])
             if drudelorentz:
                 config.materials['dispersivedtype'] = config.dtypes['complex']
@@ -163,11 +175,12 @@ class ModelBuildRun:
                 config.materials['dispersiveCdtype'] = config.dtypes['C_float_or_double']
 
             # Update estimated memory (RAM) usage
-            G.memoryusage += int(3 * config.materials['maxpoles'] *
-                                 (G.nx + 1) * (G.ny + 1) * (G.nz + 1) *
-                                 np.dtype(config.materials['dispersivedtype']).itemsize)
+            config.model_configs[G.model_num].memoryusage += int(3 *
+                                                                 config.materials['maxpoles'] *
+                                                                 (G.nx + 1) * (G.ny + 1) * (G.nz + 1) *
+                                                                 np.dtype(config.materials['dispersivedtype']).itemsize)
             G.memory_check()
-            log.info(f'\nMemory (RAM) required - updated (dispersive): ~{human_size(G.memoryusage)}\n')
+            log.info(f'\nMemory (RAM) required - updated (dispersive): ~{human_size(config.model_configs[G.model_num].memoryusage)}\n')
 
             for gb in gridbuilders:
                 gb.grid.initialise_dispersive_arrays(config.materials['dispersivedtype'])
@@ -190,103 +203,93 @@ class ModelBuildRun:
         # Check to see if numerical dispersion might be a problem
         results = dispersion_analysis(G)
         if results['error']:
-            log.warning(Fore.RED + f"\nNumerical dispersion analysis not carried \
-                        out as {results['error']}" + Style.RESET_ALL)
-        elif results['N'] < config.numdispersion['mingridsampling']:
-            raise GeneralError(f"Non-physical wave propagation: Material \
-                               '{results['material'].ID}' has wavelength sampled \
-                               by {results['N']} cells, less than required \
-                               minimum for physical wave propagation. Maximum \
-                               significant frequency estimated as {results['maxfreq']:g}Hz")
-        elif results['deltavp'] and np.abs(results['deltavp']) > config.numdispersion['maxnumericaldisp']:
-            log.warning(Fore.RED + f"\nPotentially significant numerical dispersion. \
-                        Estimated largest physical phase-velocity error is \
-                        {results['deltavp']:.2f}% in material '{results['material'].ID}' \
-                        whose wavelength sampled by {results['N']} cells. Maximum \
-                        significant frequency estimated as {results['maxfreq']:g}Hz" + Style.RESET_ALL)
+            log.warning(Fore.RED + f"\nNumerical dispersion analysis not carried out as {results['error']}" + Style.RESET_ALL)
+        elif results['N'] < config.model_configs[G.model_num].numdispersion['mingridsampling']:
+            raise GeneralError(f"Non-physical wave propagation: Material '{results['material'].ID}' has wavelength sampled by {results['N']} cells, less than required minimum for physical wave propagation. Maximum significant frequency estimated as {results['maxfreq']:g}Hz")
+        elif (results['deltavp'] and np.abs(results['deltavp']) >
+              config.model_configs[G.model_num].numdispersion['maxnumericaldisp']):
+            log.warning(Fore.RED + f"\nPotentially significant numerical dispersion. Estimated largest physical phase-velocity error is {results['deltavp']:.2f}% in material '{results['material'].ID}' whose wavelength sampled by {results['N']} cells. Maximum significant frequency estimated as {results['maxfreq']:g}Hz" + Style.RESET_ALL)
         elif results['deltavp']:
-            log.info(f"\nNumerical dispersion analysis: estimated largest physical \
-                     phase-velocity error is {results['deltavp']:.2f}% in material \
-                     '{results['material'].ID}' whose wavelength sampled by {results['N']} \
-                     cells. Maximum significant frequency estimated as {results['maxfreq']:g}Hz")
+            log.info(f"\nNumerical dispersion analysis: estimated largest physical phase-velocity error is {results['deltavp']:.2f}% in material '{results['material'].ID}' whose wavelength sampled by {results['N']} cells. Maximum significant frequency estimated as {results['maxfreq']:g}Hz")
 
     def reuse_geometry(self):
-        G = self.G
         # Reset iteration number
-        G.iteration = 0
-        config.model_configs[self.i].inputfilestr = f'\n--- Model {config.model_configs[self.i].appendmodelnumber}/{self.sim_config.model_end}, \
-                                         input file (not re-processed, i.e. geometry fixed): {self.sim_config.input_file_path}'
-        log.info(Fore.GREEN + f"{config.model_configs[self.i].inputfilestr} {'-' * (get_terminal_width() - 1 - len(config.model_configs[self.i].inputfilestr))}" + Style.RESET_ALL)
-        for grid in [G] + G.subgrids:
+        self.G.iteration = 0
+        config.model_configs[self.G.model_num].inputfilestr = f'\n--- Model {config.model_configs[self.G.model_num].appendmodelnumber}/{self.sim_config.model_end}, input file (not re-processed, i.e. geometry fixed): {self.sim_config.input_file_path}'
+        log.info(Fore.GREEN + f"{config.model_configs[self.G.model_num].inputfilestr} {'-' * (get_terminal_width() - 1 - len(config.model_configs[self.G.model_num].inputfilestr))}" + Style.RESET_ALL)
+        for grid in [self.G] + self.G.subgrids:
             grid.reset_fields()
 
     def build_scene(self):
-        G = self.G
         # API for multiple scenes / model runs
-        scene = config.model_configs[self.i].get_scene()
+        scene = config.model_configs[self.G.model_num].get_scene()
 
         # If there is no scene - process the hash commands instead
         if not scene:
             scene = Scene()
             # Parse the input file into user objects and add them to the scene
-            scene = parse_hash_commands(config.model_configs[self.i], G, scene)
+            scene = parse_hash_commands(config.model_configs[self.G.model_num], self.G, scene)
 
         # Creates the internal simulation objects.
-        scene.create_internal_objects(G)
+        scene.create_internal_objects(self.G)
 
         return scene
 
     def create_output_directory(self):
-        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[self.i].output_file_path_ext = Path(self.G.outputdirectory, config.model_configs[self.i].output_file_path_ext)
+        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).
         """
 
-        G = self.G
-
         # Write an output file in HDF5 format
-        write_hdf5_outputfile(config.model_configs[self.i].output_file_path_ext, G)
+        write_hdf5_outputfile(config.model_configs[self.G.model_num].output_file_path_ext, self.G)
 
         # Write any snapshots to file
-        if G.snapshots:
+        if self.G.snapshots:
             # Create directory and construct filename from user-supplied name
             # and model run number
-            snapshotdir = config.model_configs[self.i].snapshot_dir
+            snapshotdir = config.model_configs[self.G.model_num].snapshot_dir
             if not os.path.exists(snapshotdir):
                 os.mkdir(snapshotdir)
 
             log.info('')
-            for i, snap in enumerate(G.snapshots):
-                fn = snapshotdir / Path(config.model_configs[self.i].output_file_path.stem + '_' + snap.basefilename)
+            for i, snap in enumerate(self.G.snapshots):
+                fn = snapshotdir / Path(config.model_configs[self.G.model_num].output_file_path.stem + '_' + snap.basefilename)
                 snap.filename = fn.with_suffix('.vti')
-                pbar = tqdm(total=snap.vtkdatawritesize, leave=True, unit='byte', unit_scale=True, desc='Writing snapshot file {} of {}, {}'.format(i + 1, len(G.snapshots), os.path.split(snap.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars'])
-                snap.write_vtk_imagedata(pbar, G)
+                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,
+                            disable=not config.general['progressbars'])
+                snap.write_vtk_imagedata(pbar, self.G)
                 pbar.close()
             log.info('')
 
-    def print_resource_info(self):
-        """Print resource information on runtime and memory usage."""
+    def print_resource_info(self, tsolve):
+        """Print resource information on runtime and memory usage.
+
+        Args:
+            tsolve (float): time taken to execute solving (seconds).
+        """
 
         memGPU = ''
-        if config.cuda['gpus']:
+        if config.sim_config.general['cuda']:
             memGPU = f' host + ~{human_size(self.solver.get_memsolve())} GPU'
 
         log.info(f'\nMemory (RAM) used: ~{human_size(self.p.memory_full_info().uss)}{memGPU}')
         log.info(f'Solving time [HH:MM:SS]: {datetime.timedelta(seconds=tsolve)}')
 
     def solve(self, solver):
-        """
-        Solve using FDTD method.
+        """Solve using FDTD method.
 
         Args:
             solver (Solver): solver object.
@@ -295,25 +298,34 @@ class ModelBuildRun:
             tsolve (float): time taken to execute solving (seconds).
         """
 
-        G = self.G
-
-        if config.is_messages():
-            iterator = tqdm(range(G.iterations), desc='Running simulation, model ' + str(config.model_configs[self.i]
-                            .i + 1) + '/' + str(self.sim_config.model_end), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.general['progressbars'])
-        else:
-            iterator = range(0, G.iterations)
-
         self.create_output_directory()
-        log.info(f'Output file: {config.model_configs[self.i].output_file_path_ext}\n')
+        log.info(f'Output file: {config.model_configs[self.G.model_num].output_file_path_ext}')
+
+        # Set and check number of OpenMP threads
+        if config.sim_config.general['cpu']:
+            config.model_configs[self.G.model_num].ompthreads = set_omp_threads(config.model_configs[self.G.model_num].ompthreads)
+            log.info(f'CPU (OpenMP) threads for solving: {config.model_configs[self.G.model_num].ompthreads}')
+            if config.model_configs[self.G.model_num].ompthreads > config.sim_config.hostinfo['physicalcores']:
+                log.warning(Fore.RED + f"You have specified more threads ({config.model_configs[self.G.model_num].ompthreads}) than available physical CPU cores ({config.sim_config.hostinfo['physicalcores']}). This may lead to degraded performance." + Style.RESET_ALL)
+
+        # Print information about any GPU in use
+        if config.sim_config.general['cuda']:
+            log.info(f"GPU for solving: {config.model_configs[self.G.model_num].cuda['gpu'].deviceID} - {config.model_configs[self.G.model_num].cuda['gpu'].name}")
+
+        # Prepare iterator
+        if config.sim_config.is_messages():
+            iterator = tqdm(range(self.G.iterations), desc='\nRunning simulation, model ' + str(self.G.model_num + 1) + '/' + str(config.sim_config.model_end), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.sim_config.general['progressbars'])
+        else:
+            iterator = range(0, self.G.iterations)
 
         # Run solver
-        tsolve = self.solve(solver)
+        tsolve = solver.solve(iterator)
 
         # Write output data, i.e. field data for receivers and snapshots to file(s)
-        self.write_output_data
+        self.write_output_data()
 
         # Print resource information on runtime and memory usage
-        self.print_resource_info
+        self.print_resource_info(tsolve)
 
         return tsolve
 
@@ -326,7 +338,10 @@ class GridBuilder:
         grid = self.grid
 
         # Build the PMLS
-        pbar = tqdm(total=sum(1 for value in grid.pmlthickness.values() if value > 0), desc='Building {} Grid PML boundaries'.format(self.grid.name), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.sim_config.general['progressbars'])
+        pbar = tqdm(total=sum(1 for value in grid.pmlthickness.values() if value > 0),
+                    desc=f'Building {self.grid.name} Grid PML boundaries',
+                    ncols=get_terminal_width() - 1, file=sys.stdout,
+                    disable=not config.sim_config.general['progressbars'])
         for pml_id, thickness in grid.pmlthickness.items():
             if thickness > 0:
                 build_pml(grid, pml_id, thickness)
@@ -337,7 +352,9 @@ class GridBuilder:
         # Build the model, i.e. set the material properties (ID) for every edge
         # of every Yee cell
         log.info('')
-        pbar = tqdm(total=2, desc='Building {} Grid'.format(self.grid.name), ncols=get_terminal_width() - 1, file=sys.stdout, disable=not config.sim_config.general['progressbars'])
+        pbar = tqdm(total=2, desc=f'Building {self.grid.name} Grid',
+                    ncols=get_terminal_width() - 1, file=sys.stdout,
+                    disable=not config.sim_config.general['progressbars'])
         build_electric_components(self.grid.solid, self.grid.rigidE, self.grid.ID, self.grid)
         pbar.update()
         build_magnetic_components(self.grid.solid, self.grid.rigidH, self.grid.ID, self.grid)

gprMax/pml.py

@@ -207,7 +207,7 @@ class PML:
 
         self.CFS = G.cfs
 
-        if G.gpu is None:
+        if not config.sim_config.general['cuda']:
             self.initialise_field_arrays()
 
     def initialise_field_arrays(self):

gprMax/scene.py

@@ -16,6 +16,8 @@
 # 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 .cmds_geometry.cmds_geometry import UserObjectGeometry
 from .cmds_geometry.fractal_box_builder import FractalBoxBuilder
 from .cmds_multiple import UserObjectMulti
@@ -29,6 +31,8 @@ from .subgrids.user_objects import SubGridBase as SubGridUserBase
 from .user_inputs import create_user_input_points
 from .utilities import human_size
 
+log = logging.getLogger(__name__)
+
 
 class Scene:
     """Scene stores all of the user created objects."""
@@ -90,12 +94,11 @@ class Scene:
             try:
                 obj.create(grid, uip)
             except CmdInputError:
-                logger.exception('Error creating user input object')
+                log.exception('Error creating user input object')
 
         return self
 
     def process_singlecmds(self, G):
-
         # Check for duplicate commands and warn user if they exist
         cmds_unique = list(set(self.single_cmds))
         if len(cmds_unique) != len(self.single_cmds):
@@ -131,7 +134,8 @@ class Scene:
         self.process_cmds(self.multiple_cmds, G)
 
         # Estimate and check memory (RAM) usage
-        G.memory_check()
+        log.debug('Fix memory checks')
+        # G.memory_check()
         #snapshot_memory_check(G)
 
         # Initialise an array for volumetric material IDs (solid), boolean

gprMax/solvers.py

@@ -24,19 +24,22 @@ from .updates import CPUUpdates
 from .updates import CUDAUpdates
 
 
-def create_G():
+def create_G(model_num):
     """Create grid object according to solver.
 
+    Args:
+        model_num (int): Model number.
+
     Returns:
         G (FDTDGrid): Holds essential parameters describing the model.
     """
 
-    if config.sim_config.general['cuda']:
-        G = CUDAGrid()
+    if config.sim_config.general['cpu']:
+        G = FDTDGrid(model_num)
+    elif config.sim_config.general['cuda']:
+        G = CUDAGrid(model_num)
     elif config.sim_config.subgrid:
-        G = FDTDGrid()
-    else:
-        G = FDTDGrid()
+        G = FDTDGrid(model_num)
 
     return G
 
@@ -51,7 +54,12 @@ def create_solver(G):
         solver (Solver): solver object.
     """
 
-    if config.sim_config.gpu:
+    if config.sim_config.general['cpu']:
+        updates = CPUUpdates(G)
+        solver = Solver(updates)
+        props = updates.adapt_dispersive_config()
+        updates.set_dispersive_updates(props)
+    elif config.sim_config.general['cuda']:
         updates = CUDAUpdates(G)
         solver = Solver(updates)
     elif config.sim_config.subgrid:
@@ -60,12 +68,7 @@ def create_solver(G):
         # A large range of different functions exist to advance the time step for
         # dispersive materials. The correct function is set here based on the
         # the required numerical precision and dispersive material type.
-        props = updates.adapt_dispersive_config(config)
-        updates.set_dispersive_updates(props)
-    else:
-        updates = CPUUpdates(G)
-        solver = Solver(updates)
-        props = updates.adapt_dispersive_config(config)
+        props = updates.adapt_dispersive_config()
         updates.set_dispersive_updates(props)
 
     return solver
@@ -84,9 +87,6 @@ class Solver:
         self.updates = updates
         self.hsg = hsg
 
-    def get_G(self):
-        return self.updates.G
-
     def solve(self, iterator):
         """Time step the FDTD model.
 

gprMax/updates.py

@@ -69,7 +69,7 @@ class CPUUpdates:
         update_magnetic(self.grid.nx,
                         self.grid.ny,
                         self.grid.nz,
-                        config.hostinfo['ompthreads'],
+                        config.sim_config.hostinfo['ompthreads'],
                         self.grid.updatecoeffsH,
                         self.grid.ID,
                         self.grid.Ex,
@@ -98,11 +98,11 @@ class CPUUpdates:
     def update_electric_a(self):
         """Update electric field components."""
         # All materials are non-dispersive so do standard update.
-        if config.materials['maxpoles'] == 0:
+        if config.model_configs[self.grid.model_num].materials['maxpoles'] == 0:
             update_electric(self.grid.nx,
                             self.grid.ny,
                             self.grid.nz,
-                            config.hostinfo['ompthreads'],
+                            config.sim_config.hostinfo['ompthreads'],
                             self.grid.updatecoeffsE,
                             self.grid.ID,
                             self.grid.Ex,
@@ -118,8 +118,8 @@ class CPUUpdates:
             self.dispersive_update_a(self.grid.nx,
                                      self.grid.ny,
                                      self.grid.nz,
-                                     config.hostinfo['ompthreads'],
-                                     config.materials['maxpoles'],
+                                     config.sim_config.hostinfo['ompthreads'],
+                                     config.model_configs[self.grid.model_num].materials['maxpoles'],
                                      self.grid.updatecoeffsE,
                                      self.grid.updatecoeffsdispersive,
                                      self.grid.ID,
@@ -159,12 +159,12 @@ class CPUUpdates:
             updated after the electric field has been updated by the PML and
             source updates.
         """
-        if config.materials['maxpoles'] != 0:
+        if config.model_configs[self.grid.model_num].materials['maxpoles'] != 0:
             self.dispersive_update_b(self.grid.nx,
                                      self.grid.ny,
                                      self.grid.nz,
-                                     config.hostinfo['ompthreads'],
-                                     config.materials['maxpoles'],
+                                     config.sim_config.hostinfo['ompthreads'],
+                                     config.model_configs[self.grid.model_num].materials['maxpoles'],
                                      self.grid.updatecoeffsdispersive,
                                      self.grid.ID,
                                      self.grid.Tx,
@@ -174,26 +174,23 @@ class CPUUpdates:
                                      self.grid.Ey,
                                      self.grid.Ez)
 
-    def adapt_dispersive_config(self, config):
+    def adapt_dispersive_config(self):
         """Set properties for disperive materials.
 
-        Args:
-            config ():
-
         Returns:
             props (Props): Dispersive material properties.
         """
-        if config.materials['maxpoles'] > 1:
+        if config.model_configs[self.grid.model_num].materials['maxpoles'] > 1:
             poles = 'multi'
         else:
             poles = '1'
 
-        if config.precision == 'single':
+        if config.sim_config.general['precision'] == 'single':
             type = 'float'
         else:
             type = 'double'
 
-        if config.materials['dispersivedtype'] == config.dtypes['complex']:
+        if config.model_configs[self.grid.model_num].materials['dispersivedtype'] == config.sim_config.dtypes['complex']:
             dispersion = 'complex'
         else:
             dispersion = 'real'
@@ -276,7 +273,7 @@ class CUDAUpdates:
         """Electric and magnetic field updates - prepare kernels, and
             get kernel functions.
         """
-        if config.materials['maxpoles'] > 0:
+        if config.model_configs[self.grid.model_num].materials['maxpoles'] > 0:
             kernels_fields = SourceModule(kernels_template_fields.substitute(
                                           REAL=cudafloattype,
                                           COMPLEX=cudacomplextype,
@@ -319,7 +316,7 @@ class CUDAUpdates:
         self.copy_mat_coeffs()
 
         # Electric and magnetic field updates - dispersive materials - get kernel functions and initialise array on GPU
-        if config.materials['maxpoles'] > 0:  # If there are any dispersive materials (updates are split into two parts as they require present and updated electric field values).
+        if config.model_configs[self.grid.model_num].materials['maxpoles'] > 0:  # If there are any dispersive materials (updates are split into two parts as they require present and updated electric field values).
             self.dispersive_update_a = kernels_fields.get_function("update_electric_dispersive_A")
             self.dispersive_update_b = kernels_fields.get_function("update_electric_dispersive_B")
             self.grid.gpu_initialise_dispersive_arrays()
@@ -543,7 +540,7 @@ class CUDAUpdates:
     def update_electric_a(self):
         """Update electric field components."""
         # All materials are non-dispersive so do standard update.
-        if config.materials['maxpoles'] == 0:
+        if config.model_configs[self.grid.model_num].materials['maxpoles'] == 0:
             self.update_electric(np.int32(self.grid.nx),
                                  np.int32(self.grid.ny),
                                  np.int32(self.grid.nz),
@@ -563,7 +560,7 @@ class CUDAUpdates:
             self.dispersive_update_a(np.int32(self.grid.nx),
                                      np.int32(self.grid.ny),
                                      np.int32(self.grid.nz),
-                                     np.int32(config.materials['maxpoles']),
+                                     np.int32(config.model_configs[self.grid.model_num].materials['maxpoles']),
                                      self.grid.updatecoeffsdispersive_gpu.gpudata,
                                      self.grid.Tx_gpu,
                                      self.grid.Ty_gpu,
@@ -628,11 +625,11 @@ class CUDAUpdates:
             updated after the electric field has been updated by the PML and
             source updates.
         """
-        if config.materials['maxpoles'] != 0:
+        if config.model_configs[self.grid.model_num].materials['maxpoles'] != 0:
             self.dispersive_update_b(np.int32(self.grid.nx),
                                      np.int32(self.grid.ny),
                                      np.int32(self.grid.nz),
-                                     np.int32(config.materials['maxpoles']),
+                                     np.int32(config.model_configs[self.grid.model_num].materials['maxpoles']),
                                      self.grid.updatecoeffsdispersive_gpu.gpudata,
                                      self.grid.Tx_gpu,
                                      self.grid.Ty_gpu,

gprMax/utilities.py

@@ -340,7 +340,7 @@ def get_host_info():
     return hostinfo
 
 
-def set_openmp_threads():
+def set_omp_threads(nthreads=None):
     """Sets the number of OpenMP CPU threads for parallelised parts of code.
 
     Returns:
@@ -367,7 +367,9 @@ def set_openmp_threads():
         del os.environ['OMP_PLACES']
         del os.environ['OMP_PROC_BIND']
 
-    if os.environ.get('OMP_NUM_THREADS'):
+    if nthreads:
+        os.environ['OMP_NUM_THREADS'] = str(nthreads)
+    elif os.environ.get('OMP_NUM_THREADS'):
         nthreads = int(os.environ.get('OMP_NUM_THREADS'))
     else:
         # Set number of threads to number of physical CPU cores