Significant restructuring to handle re-using geometry for multiple model runs.

gprMax/gprMax.py

@@ -53,6 +53,7 @@ def main():
     parser.add_argument('-mpi', action='store_true', default=False, help='switch on MPI task farm')
     parser.add_argument('-benchmark', action='store_true', default=False, help='switch on benchmarking mode')
     parser.add_argument('--geometry-only', action='store_true', default=False, help='only build model and produce geometry file(s)')
+    parser.add_argument('--geometry-fixed', action='store_true', default=False, help='do not reprocess model geometry for multiple model runs')
     parser.add_argument('--write-processed', action='store_true', default=False, help='write an input file after any Python code and include commands in the original input file have been processed')
     parser.add_argument('--opt-taguchi', action='store_true', default=False, help='optimise parameters using the Taguchi optimisation method')
     args = parser.parse_args()
@@ -61,7 +62,7 @@ def main():
     inputfile = os.path.abspath(os.path.join(inputdirectory, os.path.basename(args.inputfile)))
     
     # Create a separate namespace that users can access in any Python code blocks in the input file
-    usernamespace = {'c': c, 'e0': e0, 'm0': m0, 'z0': z0, 'number_model_runs': numbermodelruns, 'inputdirectory': inputdirectory}
+    usernamespace = {'c': c, 'e0': e0, 'm0': m0, 'z0': z0, 'number_model_runs': numbermodelruns, 'input_directory': inputdirectory}
 
     # Process for Taguchi optimisation
     if args.opt_taguchi:
@@ -239,26 +240,31 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
     # Monitor memory usage
     p = psutil.Process()
     
+    # Declare variable to hold FDTDGrid class
+    global G
+    
+    # Normal model reading/building process; bypassed if geometry information to be reused
+    if not 'G' in globals():
         print('\n{}\n\nModel input file: {}\n'.format(68*'*', inputfile))
         
         # Add the current model run to namespace that can be accessed by user in any Python code blocks in input file
         usernamespace['current_model_run'] = modelrun
         print('Constants/variables available for Python scripting: {}\n'.format(usernamespace))
         
-    # Process any user input Python commands
+        # Read input file and process any Python or include commands
         processedlines = process_python_include_code(inputfile, usernamespace)
         
-    # Write a file containing the input commands after Python blocks have been processed
+        # Write a file containing the input commands after Python or include commands have been processed
         if args.write_processed:
             write_processed_file(inputfile, modelrun, numbermodelruns, processedlines)
         
-    # Check validity of command names & that essential commands are present
+        # Check validity of command names and that essential commands are present
         singlecmds, multicmds, geometry = check_cmd_names(processedlines)
 
         # Initialise an instance of the FDTDGrid class
         G = FDTDGrid()
         G.inputfilename = os.path.split(inputfile)[1]
-    G.inputdirectory = usernamespace['inputdirectory']
+        G.inputdirectory = os.path.dirname(os.path.abspath(inputfile))
 
         # Create built-in materials
         m = Material(0, 'pec', G)
@@ -274,10 +280,13 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
         process_multicmds(multicmds, G)
 
         # Initialise an array for volumetric material IDs (solid), boolean arrays for specifying materials not to be averaged (rigid),
-    # an array for cell edge IDs (ID), and arrays for the field components.
-    G.initialise_std_arrays()
+        # an array for cell edge IDs (ID)
+        G.initialise_geometry_arrays()
 
-    # Process the geometry commands in the order they were given
+        # Initialise arrays for the field components
+        G.initialise_field_arrays()
+
+        # Process geometry commands in the order they were given
         tinputprocstart = perf_counter()
         process_geometrycmds(geometry, G)
         tinputprocend = perf_counter()
@@ -298,7 +307,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
             voltagesource.create_material(G)
         
         # Initialise arrays of update coefficients to pass to update functions
-    G.initialise_std_updatecoeff_arrays()
+        G.initialise_std_update_coeff_arrays()
 
         # Initialise arrays of update coefficients and temporary values if there are any dispersive materials
         if Material.maxpoles != 0:
@@ -339,9 +348,36 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
 
         # Check to see if numerical dispersion might be a problem
         resolution = dispersion_check(G)
-    if resolution != 0 and max((G.dx, G.dy, G.dz)) > resolution:
+        if resolution and max((G.dx, G.dy, G.dz)) > resolution:
             print('\nWARNING: Potential numerical dispersion in the simulation. Check the spatial discretisation against the smallest wavelength present. Suggested resolution should be less than {:g}m'.format(resolution))
 
+    # If geometry information to be reused between model runs
+    else:        
+        # Clear arrays for field components
+        G.initialise_field_arrays()
+        
+        # Clear arrays for fields in PML
+        for pml in G.pmls:
+            pml.initialise_field_arrays()
+
+    # Adjust position of simple sources and receivers if required
+    if G.srcstepx > 0 or G.srcstepy > 0 or G.srcstepz > 0:
+        for source in itertools.chain(G.hertziandipoles, G.magneticdipoles):
+            if modelrun == 1:
+                if source.xcoord + G.srcstepx * (numbermodelruns - 1) > G.nx or source.ycoord + G.srcstepy * (numbermodelruns - 1) > G.ny or source.zcoord + G.srcstepz * (numbermodelruns - 1) > G.nz:
+                    raise GeneralError('Source(s) will be stepped to a position outside the domain.')
+            source.xcoord = source.xcoordbase + (modelrun - 1) * G.srcstepx
+            source.ycoord = source.ycoordbase + (modelrun - 1) * G.srcstepy
+            source.zcoord = source.zcoordbase + (modelrun - 1) * G.srcstepz
+    if G.rxstepx > 0 or G.rxstepy > 0 or G.rxstepz > 0:
+        for receiver in G.rxs:
+            if modelrun == 1:
+                if receiver.xcoord + G.rxstepx * (numbermodelruns - 1) > G.nx or receiver.ycoord + G.rxstepy * (numbermodelruns - 1) > G.ny or receiver.zcoord + G.rxstepz * (numbermodelruns - 1) > G.nz:
+                    raise GeneralError('Receiver(s) will be stepped to a position outside the domain.')
+            receiver.xcoord = receiver.xcoordbase + (modelrun - 1) * G.rxstepx
+            receiver.ycoord = receiver.ycoordbase + (modelrun - 1) * G.rxstepy
+            receiver.zcoord = receiver.zcoordbase + (modelrun - 1) * G.rxstepz
+    
     # Write files for any geometry views
     if not G.geometryviews and args.geometry_only:
         raise GeneralError('No geometry views found.')
@@ -352,31 +388,13 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
         tgeoend = perf_counter()
         print('\nGeometry file(s) written in [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=int(tgeoend - tgeostart))))
 
-    # Run simulation if not doing only geometry
+    # Run simulation (if not doing geometry only)
     if not args.geometry_only:
         
         # Prepare any snapshot files
         for snapshot in G.snapshots:
             snapshot.prepare_vtk_imagedata(modelrun, numbermodelruns, G)
         
-        # Adjust position of sources and receivers if required
-        if G.srcstepx > 0 or G.srcstepy > 0 or G.srcstepz > 0:
-            for source in itertools.chain(G.hertziandipoles, G.magneticdipoles, G.voltagesources, G.transmissionlines):
-                if modelrun == 1:
-                    if source.xcoord + G.srcstepx * (numbermodelruns - 1) > G.nx or source.ycoord + G.srcstepy * (numbermodelruns - 1) > G.ny or source.zcoord + G.srcstepz * (numbermodelruns - 1) > G.nz:
-                        raise GeneralError('Source(s) will be stepped to a position outside the domain.')
-                source.xcoord += (modelrun - 1) * G.srcstepx
-                source.ycoord += (modelrun - 1) * G.srcstepy
-                source.zcoord += (modelrun - 1) * G.srcstepz
-        if G.rxstepx > 0 or G.rxstepy > 0 or G.rxstepz > 0:
-            for receiver in G.rxs:
-                if modelrun == 1:
-                    if receiver.xcoord + G.rxstepx * (numbermodelruns - 1) > G.nx or receiver.ycoord + G.rxstepy * (numbermodelruns - 1) > G.ny or receiver.zcoord + G.rxstepz * (numbermodelruns - 1) > G.nz:
-                        raise GeneralError('Receiver(s) will be stepped to a position outside the domain.')
-                receiver.xcoord += (modelrun - 1) * G.rxstepx
-                receiver.ycoord += (modelrun - 1) * G.rxstepy
-                receiver.zcoord += (modelrun - 1) * G.rxstepz
-
         # Prepare output file
         inputfileparts = os.path.splitext(inputfile)
         if numbermodelruns == 1:
@@ -466,6 +484,10 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
         print('\n\nSolving took [HH:MM:SS]: {}'.format(datetime.timedelta(seconds=int(tsolveend - tsolvestart))))
         print('Peak memory (approx) used: {}'.format(human_size(p.memory_info().rss)))
 
+        # If geometry information to be reused between model runs then FDTDGrid class instance must be global so that it persists
+        if not args.geometry_fixed:
+            del G
+
         return int(tsolveend - tsolvestart)