diff --git a/gprMax/cmds_single_use.py b/gprMax/cmds_single_use.py index 2319ecdf..9b265c48 100644 --- a/gprMax/cmds_single_use.py +++ b/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): diff --git a/gprMax/config.py b/gprMax/config.py index 8a2c151a..eefadb99 100644 --- a/gprMax/config.py +++ b/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 . +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 diff --git a/gprMax/contexts.py b/gprMax/contexts.py index 72263d83..f4320349 100644 --- a/gprMax/contexts.py +++ b/gprMax/contexts.py @@ -17,6 +17,7 @@ # along with gprMax. If not, see . 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. diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py index 56321c9b..eb35b7e4 100644 --- a/gprMax/geometry_outputs.py +++ b/gprMax/geometry_outputs.py @@ -17,6 +17,7 @@ # along with gprMax. If not, see . 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): """ diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index 19e6923e..c0b439ec 100644 --- a/gprMax/gprMax.py +++ b/gprMax/gprMax.py @@ -140,4 +140,3 @@ def run_main(args): write_simulation_config(args) context = create_context() context.run() - context.print_time_report() diff --git a/gprMax/grid.py b/gprMax/grid.py index 522c9c3f..501b69d9 100644 --- a/gprMax/grid.py +++ b/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']))) diff --git a/gprMax/materials.py b/gprMax/materials.py index 8b5f9cac..fe5689f6 100644 --- a/gprMax/materials.py +++ b/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) diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py index f38b0e27..7f659aa2 100644 --- a/gprMax/model_build_run.py +++ b/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) diff --git a/gprMax/pml.py b/gprMax/pml.py index e4716ca5..38ee74f9 100644 --- a/gprMax/pml.py +++ b/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): diff --git a/gprMax/scene.py b/gprMax/scene.py index e7c5f092..2b306670 100644 --- a/gprMax/scene.py +++ b/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 . +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 diff --git a/gprMax/solvers.py b/gprMax/solvers.py index 9746f688..cfb31bb6 100644 --- a/gprMax/solvers.py +++ b/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. diff --git a/gprMax/updates.py b/gprMax/updates.py index 489ef16a..dc7ca6df 100644 --- a/gprMax/updates.py +++ b/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, diff --git a/gprMax/utilities.py b/gprMax/utilities.py index dc1b1b16..c49c5f3c 100644 --- a/gprMax/utilities.py +++ b/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