From 7e6434ede7ec4358e5aa4151f681df0c6a857d62 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 13 Feb 2020 13:23:11 +0000 Subject: [PATCH] Make single output file for main and subgrids. --- gprMax/fields_outputs.py | 125 ++++++------- gprMax/model_build_run.py | 6 +- tools/plot_Ascan.py | 362 ++++++++++++++++++++------------------ 3 files changed, 249 insertions(+), 244 deletions(-) diff --git a/gprMax/fields_outputs.py b/gprMax/fields_outputs.py index 1a8ce898..6667e3bc 100644 --- a/gprMax/fields_outputs.py +++ b/gprMax/fields_outputs.py @@ -98,90 +98,77 @@ __global__ void store_outputs(int NRX, int iteration, const int* __restrict__ rx """) -def write_hdf5_outputfiles(outputfile, G): - if G.rxs: - write_hdf5_main_grid_outputfile(outputfile, G) +def write_hdf5_outputfile(outputfile, G): + """Write an output file in HDF5 format. + + Args: + outputfile (str): Name of the output file. + G (FDTDGrid): Parameters describing a grid in a model. + """ + + # Check for any receivers in subgrids sg_rxs = [True for sg in G.subgrids if sg.rxs] + + # Create output file and write top-level meta data + if G.rxs or sg_rxs: + f = h5py.File(outputfile, 'w') + f.attrs['gprMax'] = __version__ + f.attrs['Title'] = G.title + + # Write meta data and data for main grid + if G.rxs: + write_grid(f, G) + + # Write meta data and data for any subgrids if sg_rxs: - write_hdf5_sub_grid_outputfile(outputfile, G) + for sg in G.subgrids: + grp = f.create_group('/subgrids/' + sg.name) + write_grid(grp, sg, is_subgrid=True) + + if G.rxs or sg_rxs: + log.info(f'Written output file: {outputfile.name}') -def write_hdf5_main_grid_outputfile(outputfile, G): - """Write an output file in HDF5 format. +def write_grid(basegrp, G, is_subgrid=False): + """Write grid meta data and data to HDF5 group. Args: - outputfile (str): Name of the output file. + basegrp (dict): HDF5 group. G (FDTDGrid): Parameters describing a grid in a model. + is_subgrid (bool): Is grid instance the main grid or a subgrid. """ - write_data(outputfile, G) - log.info(f'Written output file: {outputfile.name}') - - -def write_hdf5_sub_grid_outputfile(outputfile, G): - """Write an output file in HDF5 format. - - Args: - outputfile (str): Name of the output file. - G (FDTDGrid): Parameters describing a grid in a model. - """ - - stem = outputfile.stem - suffix = outputfile.suffix - parent = outputfile.parent - - for sg in G.subgrids: - # Create an outputfile for each subgrid - fp = stem + '_' + sg.name + suffix - fp = parent / Path(fp) - - f = write_data(fp, sg) - - # Write some additional meta data about the subgrid - f.attrs['is_os_sep'] = sg.is_os_sep - f.attrs['pml_separation'] = sg.pml_separation - f.attrs['subgrid_pml_thickness'] = sg.pmlthickness['x0'] - f.attrs['filter'] = sg.filter - f.attrs['ratio'] = sg.ratio - f.attrs['interpolation'] = sg.interpolation - - log.info(f'Written output file: {fp.name}') - -def write_data(outputfile, G): - """Write an output file in HDF5 format. - - Args: - outputfile (str): Name of the output file. - G (FDTDGrid): Parameters describing a grid in a model. - - Returns: - f (file object): File object. - """ - - f = h5py.File(outputfile, 'w') - f.attrs['gprMax'] = __version__ - f.attrs['Title'] = G.title - f.attrs['Iterations'] = G.iterations - f.attrs['nx_ny_nz'] = (G.nx, G.ny, G.nz) - f.attrs['dx_dy_dz'] = (G.dx, G.dy, G.dz) - f.attrs['dt'] = G.dt + # Write meta data for grid + basegrp.attrs['Iterations'] = G.iterations + basegrp.attrs['nx_ny_nz'] = (G.nx, G.ny, G.nz) + basegrp.attrs['dx_dy_dz'] = (G.dx, G.dy, G.dz) + basegrp.attrs['dt'] = G.dt nsrc = len(G.voltagesources + G.hertziandipoles + G.magneticdipoles + G.transmissionlines) - f.attrs['nsrc'] = nsrc - f.attrs['nrx'] = len(G.rxs) - f.attrs['srcsteps'] = G.srcsteps - f.attrs['rxsteps'] = G.rxsteps + basegrp.attrs['nsrc'] = nsrc + basegrp.attrs['nrx'] = len(G.rxs) + basegrp.attrs['srcsteps'] = G.srcsteps + basegrp.attrs['rxsteps'] = G.rxsteps + + if is_subgrid: + # Write additional meta data about subgrid + basegrp.attrs['is_os_sep'] = G.is_os_sep + basegrp.attrs['pml_separation'] = G.pml_separation + basegrp.attrs['subgrid_pml_thickness'] = G.pmlthickness['x0'] + basegrp.attrs['filter'] = G.filter + basegrp.attrs['ratio'] = G.ratio + basegrp.attrs['interpolation'] = G.interpolation # Create group for sources (except transmission lines); add type and positional data attributes srclist = G.voltagesources + G.hertziandipoles + G.magneticdipoles for srcindex, src in enumerate(srclist): - grp = f.create_group('/srcs/src' + str(srcindex + 1)) + grp = basegrp.create_group('srcs/src' + str(srcindex + 1)) 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; write arrays for line voltages and currents for tlindex, tl in enumerate(G.transmissionlines): - grp = f.create_group('/tls/tl' + str(tlindex + 1)) + grp = basegrp.create_group('tls/tl' + str(tlindex + 1)) grp.attrs['Position'] = (tl.xcoord * G.dx, tl.ycoord * G.dy, tl.zcoord * G.dz) grp.attrs['Resistance'] = tl.resistance grp.attrs['dl'] = tl.dl @@ -189,17 +176,15 @@ def write_data(outputfile, G): grp['Vinc'] = tl.Vinc grp['Iinc'] = tl.Iinc # Save total voltage and current - f['/tls/tl' + str(tlindex + 1) + '/Vtotal'] = tl.Vtotal - f['/tls/tl' + str(tlindex + 1) + '/Itotal'] = tl.Itotal + basegrp['tls/tl' + str(tlindex + 1) + '/Vtotal'] = tl.Vtotal + basegrp['tls/tl' + str(tlindex + 1) + '/Itotal'] = tl.Itotal # 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)) + grp = basegrp.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: - f['/rxs/rx' + str(rxindex + 1) + '/' + output] = rx.outputs[output] - - return f + basegrp['rxs/rx' + str(rxindex + 1) + '/' + output] = rx.outputs[output] diff --git a/gprMax/model_build_run.py b/gprMax/model_build_run.py index 9a9f59ed..086296b3 100644 --- a/gprMax/model_build_run.py +++ b/gprMax/model_build_run.py @@ -35,7 +35,7 @@ import gprMax.config as config from .cython.yee_cell_build import build_electric_components from .cython.yee_cell_build import build_magnetic_components from .exceptions import GeneralError -from .fields_outputs import write_hdf5_outputfiles +from .fields_outputs import write_hdf5_outputfile from .grid import dispersion_analysis from .hash_cmds_file import parse_hash_commands from .materials import Material @@ -213,8 +213,8 @@ class ModelBuildRun: to file(s). """ - # Write an output file(s) in HDF5 format - write_hdf5_outputfiles(config.get_model_config().output_file_path_ext, self.G) + # Write an output file in HDF5 format + write_hdf5_outputfile(config.get_model_config().output_file_path_ext, self.G) # Write any snapshots to file if self.G.snapshots: diff --git a/tools/plot_Ascan.py b/tools/plot_Ascan.py index 675568db..96e82d67 100644 --- a/tools/plot_Ascan.py +++ b/tools/plot_Ascan.py @@ -45,195 +45,215 @@ def mpl_plot(filename, outputs=Rx.defaultoutputs, fft=False): file = Path(filename) - # Open output file and read some attributes + # Open output file and read iterations f = h5py.File(file, 'r') - nrx = f.attrs['nrx'] - dt = f.attrs['dt'] - iterations = f.attrs['Iterations'] - time = np.linspace(0, (iterations - 1) * dt, num=iterations) + + # Paths to grid(s) to traverse for outputs + paths = ['/'] + + # Check if any subgrids and add path(s) + is_subgrids = "/subgrids" in f + if is_subgrids: + paths = paths + ['/subgrids/' + path + '/' for path in f['/subgrids'].keys()] + + # Get number of receivers in grid(s) + nrxs = [] + for path in paths: + if f[path].attrs['nrx'] > 0: + nrxs.append(f[path].attrs['nrx']) + else: + paths.remove(path) # Check there are any receivers - if nrx == 0: + if not paths: raise CmdInputError(f'No receivers found in {file}') - # Check for single output component when doing a FFT - if fft: - if not len(outputs) == 1: - raise CmdInputError('A single output must be specified when using the -fft option') + # Loop through all grids + for path in paths: + iterations = f[path].attrs['Iterations'] + nrx = f[path].attrs['nrx'] + dt = f[path].attrs['dt'] + time = np.linspace(0, (iterations - 1) * dt, num=iterations) - # New plot for each receiver - for rx in range(1, nrx + 1): - path = '/rxs/rx' + str(rx) + '/' - availableoutputs = list(f[path].keys()) + # Check for single output component when doing a FFT + if fft: + if not len(outputs) == 1: + raise CmdInputError('A single output must be specified when using the -fft option') - # If only a single output is required, create one subplot - if len(outputs) == 1: + # New plot for each receiver + for rx in range(1, nrx + 1): + rxpath = path + 'rxs/rx' + str(rx) + '/' + availableoutputs = list(f[rxpath].keys()) - # Check for polarity of output and if requested output is in file - if outputs[0][-1] == '-': - polarity = -1 - outputtext = '-' + outputs[0][0:-1] - output = outputs[0][0:-1] - else: - polarity = 1 - outputtext = outputs[0] - output = outputs[0] + # If only a single output is required, create one subplot + if len(outputs) == 1: - if output not in availableoutputs: - raise CmdInputError(f"{output} output requested to plot, but the available output for receiver 1 is {', '.join(availableoutputs)}") - - outputdata = f[path + output][:] * polarity - - # Plotting if FFT required - if fft: - # FFT - freqs, power = fft_power(outputdata, dt) - freqmaxpower = np.where(np.isclose(power, 0))[0][0] - - # Set plotting range to -60dB from maximum power or 4 times - # frequency at maximum power - try: - pltrange = np.where(power[freqmaxpower:] < -60)[0][0] + freqmaxpower + 1 - except: - pltrange = freqmaxpower * 4 - - pltrange = np.s_[0:pltrange] - - # Plot time history of output component - fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, num='rx' + str(rx), - figsize=(20, 10), facecolor='w', - edgecolor='w') - line1 = ax1.plot(time, outputdata, 'r', lw=2, label=outputtext) - ax1.set_xlabel('Time [s]') - ax1.set_ylabel(outputtext + ' field strength [V/m]') - ax1.set_xlim([0, np.amax(time)]) - ax1.grid(which='both', axis='both', linestyle='-.') - - # Plot frequency spectra - markerline, stemlines, baseline = ax2.stem(freqs[pltrange], - power[pltrange], '-.', - use_line_collection=True) - plt.setp(baseline, 'linewidth', 0) - plt.setp(stemlines, 'color', 'r') - plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r') - line2 = ax2.plot(freqs[pltrange], power[pltrange], 'r', lw=2) - ax2.set_xlabel('Frequency [Hz]') - ax2.set_ylabel('Power [dB]') - ax2.grid(which='both', axis='both', linestyle='-.') - - # Change colours and labels for magnetic field components or currents - if 'H' in outputs[0]: - plt.setp(line1, color='g') - plt.setp(line2, color='g') - plt.setp(ax1, ylabel=outputtext + ' field strength [A/m]') - plt.setp(stemlines, 'color', 'g') - plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g') - elif 'I' in outputs[0]: - plt.setp(line1, color='b') - plt.setp(line2, color='b') - plt.setp(ax1, ylabel=outputtext + ' current [A]') - plt.setp(stemlines, 'color', 'b') - plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b') - - plt.show() - - # Plotting if no FFT required - else: - fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]', - ylabel=outputtext + ' field strength [V/m]'), - num='rx' + str(rx), figsize=(20, 10), - facecolor='w', edgecolor='w') - line = ax.plot(time, outputdata, 'r', lw=2, label=outputtext) - ax.set_xlim([0, np.amax(time)]) - # ax.set_ylim([-15, 20]) - ax.grid(which='both', axis='both', linestyle='-.') - - if 'H' in output: - plt.setp(line, color='g') - plt.setp(ax, ylabel=outputtext + ', field strength [A/m]') - elif 'I' in output: - plt.setp(line, color='b') - plt.setp(ax, ylabel=outputtext + ', current [A]') - - # If multiple outputs required, create all nine subplots and populate only the specified ones - else: - fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]'), - num='rx' + str(rx), figsize=(20, 10), - facecolor='w', edgecolor='w') - if len(outputs) == 9: - gs = gridspec.GridSpec(3, 3, hspace=0.3, wspace=0.3) - else: - gs = gridspec.GridSpec(3, 2, hspace=0.3, wspace=0.3) - - for output in outputs: # Check for polarity of output and if requested output is in file - if output[-1] == 'm': + if outputs[0][-1] == '-': polarity = -1 - outputtext = '-' + output[0:-1] - output = output[0:-1] + outputtext = '-' + outputs[0][0:-1] + output = outputs[0][0:-1] else: polarity = 1 - outputtext = output + outputtext = outputs[0] + output = outputs[0] - # Check if requested output is in file if output not in availableoutputs: - raise CmdInputError(f"Output(s) requested to plot: {', '.join(outputs)}, but available output(s) for receiver {rx} in the file: {', '.join(availableoutputs)}") + raise CmdInputError(f"{output} output requested to plot, but the available output for receiver 1 is {', '.join(availableoutputs)}") - outputdata = f[path + output][:] * polarity + outputdata = f[rxpath + output][:] * polarity - if output == 'Ex': - ax = plt.subplot(gs[0, 0]) - ax.plot(time, outputdata, 'r', lw=2, label=outputtext) - ax.set_ylabel(outputtext + ', field strength [V/m]') - # ax.set_ylim([-15, 20]) - elif output == 'Ey': - ax = plt.subplot(gs[1, 0]) - ax.plot(time, outputdata, 'r', lw=2, label=outputtext) - ax.set_ylabel(outputtext + ', field strength [V/m]') - # ax.set_ylim([-15, 20]) - elif output == 'Ez': - ax = plt.subplot(gs[2, 0]) - ax.plot(time, outputdata, 'r', lw=2, label=outputtext) - ax.set_ylabel(outputtext + ', field strength [V/m]') - # ax.set_ylim([-15, 20]) - elif output == 'Hx': - ax = plt.subplot(gs[0, 1]) - ax.plot(time, outputdata, 'g', lw=2, label=outputtext) - ax.set_ylabel(outputtext + ', field strength [A/m]') - # ax.set_ylim([-0.03, 0.03]) - elif output == 'Hy': - ax = plt.subplot(gs[1, 1]) - ax.plot(time, outputdata, 'g', lw=2, label=outputtext) - ax.set_ylabel(outputtext + ', field strength [A/m]') - # ax.set_ylim([-0.03, 0.03]) - elif output == 'Hz': - ax = plt.subplot(gs[2, 1]) - ax.plot(time, outputdata, 'g', lw=2, label=outputtext) - ax.set_ylabel(outputtext + ', field strength [A/m]') - # ax.set_ylim([-0.03, 0.03]) - elif output == 'Ix': - ax = plt.subplot(gs[0, 2]) - ax.plot(time, outputdata, 'b', lw=2, label=outputtext) - ax.set_ylabel(outputtext + ', current [A]') - elif output == 'Iy': - ax = plt.subplot(gs[1, 2]) - ax.plot(time, outputdata, 'b', lw=2, label=outputtext) - ax.set_ylabel(outputtext + ', current [A]') - elif output == 'Iz': - ax = plt.subplot(gs[2, 2]) - ax.plot(time, outputdata, 'b', lw=2, label=outputtext) - ax.set_ylabel(outputtext + ', current [A]') - for ax in fig.axes: - ax.set_xlim([0, np.amax(time)]) - ax.grid(which='both', axis='both', linestyle='-.') + # Plotting if FFT required + if fft: + # FFT + freqs, power = fft_power(outputdata, dt) + freqmaxpower = np.where(np.isclose(power, 0))[0][0] - # Save a PDF/PNG of the figure - savename = file.stem + '_rx' + str(rx) - savename = file.parent / savename - # fig.savefig(savename.with_suffix('.pdf'), dpi=None, format='pdf', - # bbox_inches='tight', pad_inches=0.1) - # fig.savefig(savename.with_suffix('.png'), dpi=150, format='png', - # bbox_inches='tight', pad_inches=0.1) + # Set plotting range to -60dB from maximum power or 4 times + # frequency at maximum power + try: + pltrange = np.where(power[freqmaxpower:] < -60)[0][0] + freqmaxpower + 1 + except: + pltrange = freqmaxpower * 4 + + pltrange = np.s_[0:pltrange] + + # Plot time history of output component + fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, + num=rxpath + ' - ' + f[rxpath].attrs['Name'], + figsize=(20, 10), facecolor='w', + edgecolor='w') + line1 = ax1.plot(time, outputdata, 'r', lw=2, label=outputtext) + ax1.set_xlabel('Time [s]') + ax1.set_ylabel(outputtext + ' field strength [V/m]') + ax1.set_xlim([0, np.amax(time)]) + ax1.grid(which='both', axis='both', linestyle='-.') + + # Plot frequency spectra + markerline, stemlines, baseline = ax2.stem(freqs[pltrange], + power[pltrange], '-.', + use_line_collection=True) + plt.setp(baseline, 'linewidth', 0) + plt.setp(stemlines, 'color', 'r') + plt.setp(markerline, 'markerfacecolor', 'r', 'markeredgecolor', 'r') + line2 = ax2.plot(freqs[pltrange], power[pltrange], 'r', lw=2) + ax2.set_xlabel('Frequency [Hz]') + ax2.set_ylabel('Power [dB]') + ax2.grid(which='both', axis='both', linestyle='-.') + + # Change colours and labels for magnetic field components or currents + if 'H' in outputs[0]: + plt.setp(line1, color='g') + plt.setp(line2, color='g') + plt.setp(ax1, ylabel=outputtext + ' field strength [A/m]') + plt.setp(stemlines, 'color', 'g') + plt.setp(markerline, 'markerfacecolor', 'g', 'markeredgecolor', 'g') + elif 'I' in outputs[0]: + plt.setp(line1, color='b') + plt.setp(line2, color='b') + plt.setp(ax1, ylabel=outputtext + ' current [A]') + plt.setp(stemlines, 'color', 'b') + plt.setp(markerline, 'markerfacecolor', 'b', 'markeredgecolor', 'b') + + plt.show() + + # Plotting if no FFT required + else: + fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]', + ylabel=outputtext + ' field strength [V/m]'), + num=rxpath + ' - ' + f[rxpath].attrs['Name'], + figsize=(20, 10), facecolor='w', edgecolor='w') + line = ax.plot(time, outputdata, 'r', lw=2, label=outputtext) + ax.set_xlim([0, np.amax(time)]) + # ax.set_ylim([-15, 20]) + ax.grid(which='both', axis='both', linestyle='-.') + + if 'H' in output: + plt.setp(line, color='g') + plt.setp(ax, ylabel=outputtext + ', field strength [A/m]') + elif 'I' in output: + plt.setp(line, color='b') + plt.setp(ax, ylabel=outputtext + ', current [A]') + + # If multiple outputs required, create all nine subplots and populate only the specified ones + else: + fig, ax = plt.subplots(subplot_kw=dict(xlabel='Time [s]'), + num=rxpath + ' - ' + f[rxpath].attrs['Name'], + figsize=(20, 10), facecolor='w', edgecolor='w') + if len(outputs) == 9: + gs = gridspec.GridSpec(3, 3, hspace=0.3, wspace=0.3) + else: + gs = gridspec.GridSpec(3, 2, hspace=0.3, wspace=0.3) + + for output in outputs: + # Check for polarity of output and if requested output is in file + if output[-1] == 'm': + polarity = -1 + outputtext = '-' + output[0:-1] + output = output[0:-1] + else: + polarity = 1 + outputtext = output + + # Check if requested output is in file + if output not in availableoutputs: + raise CmdInputError(f"Output(s) requested to plot: {', '.join(outputs)}, but available output(s) for receiver {rx} in the file: {', '.join(availableoutputs)}") + + outputdata = f[rxpath + output][:] * polarity + + if output == 'Ex': + ax = plt.subplot(gs[0, 0]) + ax.plot(time, outputdata, 'r', lw=2, label=outputtext) + ax.set_ylabel(outputtext + ', field strength [V/m]') + # ax.set_ylim([-15, 20]) + elif output == 'Ey': + ax = plt.subplot(gs[1, 0]) + ax.plot(time, outputdata, 'r', lw=2, label=outputtext) + ax.set_ylabel(outputtext + ', field strength [V/m]') + # ax.set_ylim([-15, 20]) + elif output == 'Ez': + ax = plt.subplot(gs[2, 0]) + ax.plot(time, outputdata, 'r', lw=2, label=outputtext) + ax.set_ylabel(outputtext + ', field strength [V/m]') + # ax.set_ylim([-15, 20]) + elif output == 'Hx': + ax = plt.subplot(gs[0, 1]) + ax.plot(time, outputdata, 'g', lw=2, label=outputtext) + ax.set_ylabel(outputtext + ', field strength [A/m]') + # ax.set_ylim([-0.03, 0.03]) + elif output == 'Hy': + ax = plt.subplot(gs[1, 1]) + ax.plot(time, outputdata, 'g', lw=2, label=outputtext) + ax.set_ylabel(outputtext + ', field strength [A/m]') + # ax.set_ylim([-0.03, 0.03]) + elif output == 'Hz': + ax = plt.subplot(gs[2, 1]) + ax.plot(time, outputdata, 'g', lw=2, label=outputtext) + ax.set_ylabel(outputtext + ', field strength [A/m]') + # ax.set_ylim([-0.03, 0.03]) + elif output == 'Ix': + ax = plt.subplot(gs[0, 2]) + ax.plot(time, outputdata, 'b', lw=2, label=outputtext) + ax.set_ylabel(outputtext + ', current [A]') + elif output == 'Iy': + ax = plt.subplot(gs[1, 2]) + ax.plot(time, outputdata, 'b', lw=2, label=outputtext) + ax.set_ylabel(outputtext + ', current [A]') + elif output == 'Iz': + ax = plt.subplot(gs[2, 2]) + ax.plot(time, outputdata, 'b', lw=2, label=outputtext) + ax.set_ylabel(outputtext + ', current [A]') + for ax in fig.axes: + ax.set_xlim([0, np.amax(time)]) + ax.grid(which='both', axis='both', linestyle='-.') + + # Save a PDF/PNG of the figure + savename = file.stem + '_rx' + str(rx) + savename = file.parent / savename + # fig.savefig(savename.with_suffix('.pdf'), dpi=None, format='pdf', + # bbox_inches='tight', pad_inches=0.1) + # fig.savefig(savename.with_suffix('.png'), dpi=150, format='png', + # bbox_inches='tight', pad_inches=0.1) f.close()