working up to solver.solve()

gprMax/cmds_single_use.py

@@ -19,15 +19,17 @@
 import inspect
 import os
 import sys
+import decimal as d
 
 from .config import c
 from .config import dtypes
+from .config import hostinfo
 from .exceptions import CmdInputError
 from .waveforms import Waveform
+from .utilities import round_value
 
 import numpy as np
 from scipy import interpolate
-
 from colorama import init
 from colorama import Fore
 from colorama import Style
@@ -107,16 +109,14 @@ class DomainSingle(UserObjectSingle):
             G.mode = '3D'
 
         # Round down time step to nearest float with precision one less than hardware maximum. Avoids inadvertently exceeding the CFL due to binary representation of floating point number.
-        G.round_time_step()
+        # Round down time step to nearest float with precision one less than hardware maximum.
+        # Avoids inadvertently exceeding the CFL due to binary representation of floating point number.
+        G.dt = round_value(G.dt, decimalplaces=d.getcontext().prec - 1)
 
         if G.messages:
             print('Mode: {}'.format(G.mode))
             print('Time step (at CFL limit): {:g} secs'.format(G.dt))
 
-        # do threads here
-        from .utilities import get_host_info
-        hostinfo = get_host_info()
-
         # Number of threads (OpenMP) to use
         if sys.platform == 'darwin':
             os.environ['OMP_WAIT_POLICY'] = 'ACTIVE'  # Should waiting threads consume CPU power (can drastically effect performance)
@@ -126,7 +126,7 @@ class DomainSingle(UserObjectSingle):
         # os.environ['OMP_DISPLAY_ENV'] = 'TRUE' # Prints OMP version and environment variables (useful for debug)
 
         # Catch bug with Windows Subsystem for Linux (https://github.com/Microsoft/BashOnWindows/issues/785)
-        if 'Microsoft' in G.hostinfo['osversion']:
+        if 'Microsoft' in hostinfo['osversion']:
             os.environ['KMP_AFFINITY'] = 'disabled'
             del os.environ['OMP_PLACES']
             del os.environ['OMP_PROC_BIND']
@@ -140,7 +140,7 @@ class DomainSingle(UserObjectSingle):
 
         if G.messages:
             print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
-        if G.nthreads > G.hostinfo['physicalcores']:
+        if G.nthreads > hostinfo['physicalcores']:
             print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, hostinfo['physicalcores']) + Style.RESET_ALL)
 
 
@@ -180,20 +180,21 @@ class Discretisation(UserObjectSingle):
     def create(self, G, uip):
 
         try:
+            G.dl = np.array(self.kwargs['p1'])
             G.dx, G.dy, G.dz = self.kwargs['p1']
 
         except KeyError:
             raise CmdInputError('Discretisation requires a point')
 
-        if G.dx <= 0:
+        if G.dl[0] <= 0:
             raise CmdInputError('Discretisation requires the x-direction spatial step to be greater than zero')
-        if G.dy <= 0:
+        if G.dl[1] <= 0:
             raise CmdInputError(' Discretisation requires the y-direction spatial step to be greater than zero')
-        if G.dz <= 0:
+        if G.dl[2] <= 0:
             raise CmdInputError('Discretisation requires the z-direction spatial step to be greater than zero')
 
         if G.messages:
-            print('Spatial discretisation: {:g} x {:g} x {:g}m'.format(G.dx, G.dy, G.dz))
+            print('Spatial discretisation: {:g} x {:g} x {:g}m'.format(*G.dl))
 
 
 class TimeWindow(UserObjectSingle):
@@ -307,8 +308,6 @@ class NumThreads(UserObjectSingle):
 
     def create(self, G, uip):
         # Get information about host machine
-        from .utilities import get_host_info
-        hostinfo = get_host_info()
 
         try:
             n = self.kwargs['n']
@@ -323,7 +322,7 @@ class NumThreads(UserObjectSingle):
 
         if G.messages:
             print('Number of CPU (OpenMP) threads: {}'.format(G.nthreads))
-        if G.nthreads > G.hostinfo['physicalcores']:
+        if G.nthreads > hostinfo['physicalcores']:
             print(Fore.RED + 'WARNING: You have specified more threads ({}) than available physical CPU cores ({}). This may lead to degraded performance.'.format(G.nthreads, hostinfo['physicalcores']) + Style.RESET_ALL)
 
         # Print information about any GPU in use

gprMax/config.py

@@ -46,7 +46,7 @@ z0 = np.sqrt(m0 / e0)
 #   mode: 2D TMx, 2D TMy, 2D TMz, or 3D
 #   cpu, cuda, opencl: solver type
 general = {'inputfilepath': '', 'outputfilepath': '', 'messages': True,
-           'progressbars': True, 'mode': None, 'cpu': True, 'cuda': False, 'opencl': False}
+           'progressbars': True, 'mode': '3D', 'cpu': True, 'cuda': False, 'opencl': False}
 
 # Store information about host machine
 hostinfo = get_host_info()
@@ -141,6 +141,7 @@ class SimulationConfig:
         self.general = {}
         self.general['messages'] = True
         self.geometry_fixed = args.geometry_fixed
+        self.geometry_only = args.geometry_only
         self.scenes = args.scenes
 
         self.set_input_file_path()

gprMax/grid.py

@@ -65,14 +65,38 @@ class Grid(object):
     def get(self, i, j, k):
         return self.grid[i, j, k]
 
-    def within_bounds(self, **kwargs):
+    def within_bounds(self, p):
-        for co, val in kwargs.items():
+        if p[0] < 0 or p[0] > self.nx:
-            if val < 0 or val > getattr(self, 'n' + co):
+            raise ValueError('x')
-                raise ValueError(co)
+        if p[1] < 0 or p[1] > self.ny:
+            raise ValueError('y')
+        if p[2] < 0 or p[2] > self.nz:
+            raise ValueError('z')
 
-    def calculate_coord(self, coord, val):
+    def discretise_point(self, p):
-        co = round_value(float(val) / getattr(self, 'd' + coord))
+        x = round_value(float(p[0]) / self.dx)
-        return co
+        y = round_value(float(p[1]) / self.dy)
+        z = round_value(float(p[2]) / self.dz)
+        return (x, y, z)
+
+    def round_to_grid(self, p):
+        p = self.discretise_point(p)
+        p_r = (p[0] * self.dx,
+               p[1] * self.dy,
+               p[2] * self.dz)
+        return p_r
+
+
+    def within_pml(self, p):
+        if (p[0] < self.pmlthickness['x0'] or
+            p[0] > self.nx - self.pmlthickness['xmax'] or
+            p[1] < self.pmlthickness['y0'] or
+            p[1] > self.ny - self.pmlthickness['ymax'] or
+            p[2] < self.pmlthickness['z0'] or
+            p[2] > self.nz - self.pmlthickness['zmax']):
+            return True
+        else:
+            return False
 
 
 class FDTDGrid(Grid):
@@ -119,6 +143,7 @@ class FDTDGrid(Grid):
         self.srcsteps = [0, 0, 0]
         self.rxsteps = [0, 0, 0]
         self.snapshots = []
+        self.subgrids = []
 
     def initialise_geometry_arrays(self):
         """
@@ -143,6 +168,13 @@ class FDTDGrid(Grid):
         self.Hy = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.dtypes['float_or_double'])
         self.Hz = np.zeros((self.nx + 1, self.ny + 1, self.nz + 1), dtype=config.dtypes['float_or_double'])
 
+
+    def initialise_grids(self):
+        """Function to call the initialisation of all grids."""
+        for g in [self] + self.subgrids:
+            g.initialise_geometry_arrays()
+            g.initialise_field_arrays()
+
     def initialise_std_update_coeff_arrays(self):
         """Initialise arrays for storing update coefficients."""
         self.updatecoeffsE = np.zeros((len(self.materials), 5), dtype=config.dtypes['float_or_double'])
@@ -362,9 +394,9 @@ def dispersion_analysis(G):
         minwavelength = minvelocity / results['maxfreq']
 
         # Maximum spatial step
-        if '3D' in config.mode:
+        if '3D' in config.general['mode']:
             delta = max(G.dx, G.dy, G.dz)
-        elif '2D' in config.mode:
+        elif '2D' in config.general['mode']:
             if G.nx == 1:
                 delta = max(G.dy, G.dz)
             elif G.ny == 1:

gprMax/model_build_run.py

@@ -110,6 +110,8 @@ class ModelBuildRun:
         else:
             self.reuse_geometry()
 
+        G = self.G
+
         # 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):
@@ -129,7 +131,7 @@ class ModelBuildRun:
                 receiver.zcoord = receiver.zcoordorigin + (currentmodelrun - 1) * G.rxsteps[2]
 
         # Write files for any geometry views and geometry object outputs
-        if not (G.geometryviews or G.geometryobjectswrite) and args.geometry_only and config.general['messages']:
+        if not (G.geometryviews or G.geometryobjectswrite) and self.sim_config.geometry_only and config.general['messages']:
             print(Fore.RED + '\nWARNING: No geometry views or geometry objects to output found.' + Style.RESET_ALL)
         if config.general['messages']: print()
         for i, geometryview in enumerate(G.geometryviews):
@@ -143,7 +145,7 @@ class ModelBuildRun:
             pbar.close()
 
         # If only writing geometry information
-        if args.geometry_only:
+        if self.sim_config.geometry_only:
             tsolve = 0
 
     def build_geometry(self):

gprMax/pml.py

@@ -445,7 +445,7 @@ def pml_information(G):
         pmlinfo = pmlinfo[:-2] + ' cells'
     return 'PML: formulation: {}, order: {}, thickness: {}'.format(G.pmlformulation, len(G.cfs), pmlinfo)
 
-def build_pml(G, key, thickness):
+def build_pml(G, key, value):
     """This function builds instances of the PML and calculates the initial
         parameters and coefficients including setting profile
         (based on underlying material er and mr from solid array).

gprMax/scene.py

@@ -5,7 +5,7 @@ from .cmds_single_use import DomainSingle
 from .cmds_single_use import Discretisation
 from .cmds_single_use import TimeWindow
 from .cmds_multiple import UserObjectMulti
-from .subgrids.user_objects import SubGridBase
+from .subgrids.user_objects import SubGridBase as SubGridUserBase
 from .cmds_geometry.cmds_geometry import UserObjectGeometry
 from .exceptions import CmdInputError
 from .cmds_geometry.fractal_box_builder import FractalBoxBuilder

gprMax/user_inputs.py

@@ -66,6 +66,7 @@ class MainGridUserInput(UserInput):
         super().__init__(grid)
 
     def check_point(self, p, cmd_str, name=''):
+        """Discretise point and check its within the domain"""
         p = self.discretise_point(p)
         self.point_within_bounds(p, cmd_str, name)
         return p