Receiver instances now store time history field component values, instead of writing to file on each iteration. An output file is written at the end of the simulation.

gprMax/gprMax.py

@@ -33,15 +33,16 @@ from ._version import __version__
 from .constants import c, e0, m0, z0
 from .exceptions import GeneralError
 from .fields_update import update_electric, update_magnetic, update_electric_dispersive_multipole_A, update_electric_dispersive_multipole_B, update_electric_dispersive_1pole_A, update_electric_dispersive_1pole_B
-from .grid import FDTDGrid, dispersion_check
+from .grid import FDTDGrid, dispersion_check, Ix, Iy, Iz
 from .input_cmds_geometry import process_geometrycmds
 from .input_cmds_file import process_python_include_code, write_processed_file, check_cmd_names
 from .input_cmds_multiuse import process_multicmds
 from .input_cmds_singleuse import process_singlecmds
 from .materials import Material
-from .writer_hdf5 import prepare_hdf5, write_hdf5
 from .pml import build_pmls, update_electric_pml, update_magnetic_pml
+from .receivers import store_outputs
 from .utilities import update_progress, logo, human_size
+from .writer_hdf5 import write_hdf5
 from .yee_cell_build import build_electric_components, build_magnetic_components
 
 
@@ -303,10 +304,10 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
         G.inputdirectory = os.path.dirname(os.path.abspath(inputfile))
 
         # Create built-in materials
-        m = Material(0, 'pec', G)
+        m = Material(0, 'pec')
         m.average = False
         G.materials.append(m)
-        m = Material(1, 'free_space', G)
+        m = Material(1, 'free_space')
         G.materials.append(m)
 
         # Process parameters for commands that can only occur once in the model
@@ -440,11 +441,10 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
             outputfile = inputfileparts[0] + str(modelrun) + '.out'
         sys.stdout.write('\nOutput to file: {}\n'.format(outputfile))
         sys.stdout.flush()
-        f = prepare_hdf5(outputfile, G)
 
-        ##################################
+        ####################################
-        #   Main FDTD calculation loop   #
+        #  Start - Main FDTD calculations  #
-        ##################################
+        ####################################
         tsolvestart = perf_counter()
         # Absolute time
         abstime = 0
@@ -453,8 +453,8 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
             if timestep == 0:
                 tstepstart = perf_counter()
 
-            # Write field outputs to file
+            # Store field component values for every receiver and transmission line
-            write_hdf5(f, timestep, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G)
+            store_outputs(timestep, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G)
 
             # Write any snapshots to file
             for snapshot in G.snapshots:
@@ -514,13 +514,17 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
             elif timestep > 1:
                 update_progress((timestep + 1) / G.iterations)
 
-        # Close output file
+        # Write an output file in HDF5 format
-        f.close()
+        write_hdf5(outputfile, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G)
 
         tsolveend = perf_counter()
         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)))
 
+        ##################################
+        #  End - Main FDTD calculations  #
+        ##################################
+
         # 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

gprMax/receivers.py

@@ -16,18 +16,55 @@
 # You should have received a copy of the GNU General Public License
 # along with gprMax.  If not, see <http://www.gnu.org/licenses/>.
 
+from gprMax.grid import Ix, Iy, Iz
+
+
 class Rx(object):
-    """Receiever output points."""
+    """Receiver output points."""
     
     availableoutputs = ['Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz', 'Ix', 'Iy', 'Iz']
     
     def __init__(self):
         
         self.ID = None
-        self.outputs = []
+        self.outputs = dict()
         self.xcoord = None
         self.ycoord = None
         self.zcoord = None
         self.xcoordbase = None
         self.ycoordbase = None
         self.zcoordbase = None
+
+
+def store_outputs(timestep, Ex, Ey, Ez, Hx, Hy, Hz, G):
+    """Stores field component values for every receiver and transmission line.
+
+    Args:
+        timestep (int): Current iteration number.
+        Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic field values.
+        G (class): Grid class instance - holds essential parameters describing the model.
+    """
+
+    for rx in G.rxs:
+        if 'Ex' in rx.outputs:
+            rx.outputs['Ex'][timestep] = Ex[rx.xcoord, rx.ycoord, rx.zcoord]
+        if 'Ey' in rx.outputs:
+            rx.outputs['Ey'][timestep] = Ey[rx.xcoord, rx.ycoord, rx.zcoord]
+        if 'Ez' in rx.outputs:
+            rx.outputs['Ez'][timestep] = Ez[rx.xcoord, rx.ycoord, rx.zcoord]
+        if 'Hx' in rx.outputs:
+            rx.outputs['Hx'][timestep] = Hx[rx.xcoord, rx.ycoord, rx.zcoord]
+        if 'Hy' in rx.outputs:
+            rx.outputs['Hy'][timestep] = Hy[rx.xcoord, rx.ycoord, rx.zcoord]
+        if 'Hz' in rx.outputs:
+            rx.outputs['Hz'][timestep] = Hz[rx.xcoord, rx.ycoord, rx.zcoord]
+        if 'Ix' in rx.outputs:
+            rx.outputs['Ix'][timestep] = Ix(rx.xcoord, rx.ycoord, rx.zcoord, Hy, Hz, G)
+        if 'Iy' in rx.outputs:
+            rx.outputs['Iy'][timestep] = Iy(rx.xcoord, rx.ycoord, rx.zcoord, Hx, Hz, G)
+        if 'Iz' in rx.outputs:
+            rx.outputs['Iz'][timestep] = Iz(rx.xcoord, rx.ycoord, rx.zcoord, Hx, Hy, G)
+
+    for tlindex, tl in enumerate(G.transmissionlines):
+        tl.Vtotal[timestep] = tl.voltage[tl.antpos]
+        tl.Itotal[timestep] = tl.current[tl.antpos]

gprMax/sources.py

@@ -226,6 +226,8 @@ class TransmissionLine(Source):
         self.current = np.zeros(self.nl, dtype=floattype)
         self.Vinc = np.zeros(G.iterations, dtype=floattype)
         self.Iinc = np.zeros(G.iterations, dtype=floattype)
+        self.Vtotal = np.zeros(G.iterations, dtype=floattype)
+        self.Itotal = np.zeros(G.iterations, dtype=floattype)
     
     def calculate_incident_V_I(self, G):
         """Calculates the incident voltage and current with a long length transmission line not connected to the main grid from: http://dx.doi.org/10.1002/mop.10415

gprMax/writer_hdf5.py

@@ -20,18 +20,15 @@ import h5py
 
 from gprMax._version import __version__
 from gprMax.constants import floattype
-from gprMax.grid import Ix, Iy, Iz
 
 
-def prepare_hdf5(outputfile, G):
+def write_hdf5(outputfile, Ex, Ey, Ez, Hx, Hy, Hz, G):
-    """Prepares an output file in HDF5 format for writing.
+    """Write an output file in HDF5 format.
 
     Args:
         outputfile (str): Name of the output file.
+        Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic field values.
         G (class): Grid class instance - holds essential parameters describing the model.
-
-    Returns:
-        f (file object): File object for the file to be written to.
     """
 
     f = h5py.File(outputfile, 'w')
@@ -55,7 +52,7 @@ def prepare_hdf5(outputfile, G):
         grp.attrs['Type'] = type(src).__name__
         grp.attrs['Position'] = (src.xcoord * G.dx, src.ycoord * G.dy, src.zcoord * G.dz)
 
-    # Create group for transmission lines; add positional data, line resistance and line discretisation attributes; initialise arrays for line voltages and currents
+    # Create group for transmission lines; add positional data, line resistance and line discretisation attributes; write arrays for line voltages and currents
     if G.transmissionlines:
         for tlindex, tl in enumerate(G.transmissionlines):
             grp = f.create_group('/tls/tl' + str(tlindex + 1))
@@ -65,53 +62,34 @@ def prepare_hdf5(outputfile, G):
             # Save incident voltage and current
             grp['Vinc'] = tl.Vinc
             grp['Iinc'] = tl.Iinc
-            grp.create_dataset('Vtotal', (G.iterations, ), dtype=floattype)
+            # Save total voltage and current
-            grp.create_dataset('Itotal', (G.iterations, ), dtype=floattype)
+            f['/tls/tl' + str(tlindex + 1) + '/Vtotal'] = tl.Vtotal
+            f['/tls/tl' + str(tlindex + 1) + '/Itotal'] = tl.Itotal
 
-    # Create group and add positional data and initialise field component arrays for receivers
+
+    # Create group, add positional data and write field component arrays for receivers
     for rxindex, rx in enumerate(G.rxs):
         grp = f.create_group('/rxs/rx' + str(rxindex + 1))
         if rx.ID:
             grp.attrs['Name'] = rx.ID
         grp.attrs['Position'] = (rx.xcoord * G.dx, rx.ycoord * G.dy, rx.zcoord * G.dz)
-        for output in rx.outputs:
-            grp.create_dataset(output, (G.iterations, ), dtype=floattype)
 
-    return f
-
-
-def write_hdf5(f, timestep, Ex, Ey, Ez, Hx, Hy, Hz, G):
-    """Writes field component values to an output file in HDF5 format.
-
-    Args:
-        f (file object): File object for the file to be written to.
-        timestep (int): Current iteration number.
-        Ex, Ey, Ez, Hx, Hy, Hz (memory view): Current electric and magnetic field values.
-        G (class): Grid class instance - holds essential parameters describing the model.
-    """
-
-    # For each rx, write field component values at current timestep
-    for rxindex, rx in enumerate(G.rxs):
         if 'Ex' in rx.outputs:
-            f['/rxs/rx' + str(rxindex + 1) + '/Ex'][timestep] = Ex[rx.xcoord, rx.ycoord, rx.zcoord]
+            f['/rxs/rx' + str(rxindex + 1) + '/Ex'] = rx.outputs['Ex']
         if 'Ey' in rx.outputs:
-            f['/rxs/rx' + str(rxindex + 1) + '/Ey'][timestep] = Ey[rx.xcoord, rx.ycoord, rx.zcoord]
+            f['/rxs/rx' + str(rxindex + 1) + '/Ey'] = rx.outputs['Ey']
         if 'Ez' in rx.outputs:
-            f['/rxs/rx' + str(rxindex + 1) + '/Ez'][timestep] = Ez[rx.xcoord, rx.ycoord, rx.zcoord]
+            f['/rxs/rx' + str(rxindex + 1) + '/Ez'] = rx.outputs['Ez']
         if 'Hx' in rx.outputs:
-            f['/rxs/rx' + str(rxindex + 1) + '/Hx'][timestep] = Hx[rx.xcoord, rx.ycoord, rx.zcoord]
+            f['/rxs/rx' + str(rxindex + 1) + '/Hx'] = rx.outputs['Hx']
         if 'Hy' in rx.outputs:
-            f['/rxs/rx' + str(rxindex + 1) + '/Hy'][timestep] = Hy[rx.xcoord, rx.ycoord, rx.zcoord]
+            f['/rxs/rx' + str(rxindex + 1) + '/Hy'] = rx.outputs['Hy']
         if 'Hz' in rx.outputs:
-            f['/rxs/rx' + str(rxindex + 1) + '/Hz'][timestep] = Hz[rx.xcoord, rx.ycoord, rx.zcoord]
+            f['/rxs/rx' + str(rxindex + 1) + '/Hz'] = rx.outputs['Hz']
         if 'Ix' in rx.outputs:
-            f['/rxs/rx' + str(rxindex + 1) + '/Ix'][timestep] = Ix(rx.xcoord, rx.ycoord, rx.zcoord, G.Hy, G.Hz, G)
+            f['/rxs/rx' + str(rxindex + 1) + '/Ix'] = rx.outputs['Ix']
         if 'Iy' in rx.outputs:
-            f['/rxs/rx' + str(rxindex + 1) + '/Iy'][timestep] = Iy(rx.xcoord, rx.ycoord, rx.zcoord, G.Hx, G.Hz, G)
+            f['/rxs/rx' + str(rxindex + 1) + '/Iy'] = rx.outputs['Iy']
         if 'Iz' in rx.outputs:
-            f['/rxs/rx' + str(rxindex + 1) + '/Iz'][timestep] = Iz(rx.xcoord, rx.ycoord, rx.zcoord, G.Hx, G.Hy, G)
+            f['/rxs/rx' + str(rxindex + 1) + '/Iz'] = rx.outputs['Iz']
 
-    if G.transmissionlines:
-        for tlindex, tl in enumerate(G.transmissionlines):
-            f['/tls/tl' + str(tlindex + 1) + '/Vtotal'][timestep] = tl.voltage[tl.antpos]
-            f['/tls/tl' + str(tlindex + 1) + '/Itotal'][timestep] = tl.current[tl.antpos]