From 31e8fde42ca248952d39991b638afd10320ef573 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 30 Jun 2016 14:11:18 +0100 Subject: [PATCH] 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 | 28 ++++++++++++--------- gprMax/receivers.py | 43 +++++++++++++++++++++++++++++--- gprMax/sources.py | 2 ++ gprMax/writer_hdf5.py | 58 ++++++++++++++----------------------------- 4 files changed, 76 insertions(+), 55 deletions(-) diff --git a/gprMax/gprMax.py b/gprMax/gprMax.py index 17987f1d..354048e9 100644 --- a/gprMax/gprMax.py +++ b/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 - write_hdf5(f, timestep, G.Ex, G.Ey, G.Ez, G.Hx, G.Hy, G.Hz, G) + # Store field component values for every receiver and transmission line + 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 - f.close() + # Write an output file in HDF5 format + 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 diff --git a/gprMax/receivers.py b/gprMax/receivers.py index 5ac050fd..b992f856 100644 --- a/gprMax/receivers.py +++ b/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 . +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 \ No newline at end of file + 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] \ No newline at end of file diff --git a/gprMax/sources.py b/gprMax/sources.py index c75d6101..209c5e15 100644 --- a/gprMax/sources.py +++ b/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 diff --git a/gprMax/writer_hdf5.py b/gprMax/writer_hdf5.py index 63516f3d..7f82c593 100644 --- a/gprMax/writer_hdf5.py +++ b/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): - """Prepares an output file in HDF5 format for writing. +def write_hdf5(outputfile, Ex, Ey, Ez, Hx, Hy, Hz, G): + """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) - grp.create_dataset('Itotal', (G.iterations, ), dtype=floattype) + # Save total voltage and current + 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]