From 8c845b6151e87adcf705a4fc91062a3d650b69f9 Mon Sep 17 00:00:00 2001 From: AngelFP Date: Wed, 18 Aug 2021 13:43:41 +0200 Subject: [PATCH 1/4] Implement derived `OpenPMDDataReader` class --- visualpic/data_reading/openpmd_data_reader.py | 245 ++++++++++++++++++ 1 file changed, 245 insertions(+) create mode 100644 visualpic/data_reading/openpmd_data_reader.py diff --git a/visualpic/data_reading/openpmd_data_reader.py b/visualpic/data_reading/openpmd_data_reader.py new file mode 100644 index 0000000..9e47d8f --- /dev/null +++ b/visualpic/data_reading/openpmd_data_reader.py @@ -0,0 +1,245 @@ +import h5py +from openpmd_viewer.openpmd_timeseries.data_reader import DataReader +from openpmd_viewer.openpmd_timeseries.data_reader.h5py_reader import ( + field_reader as fr) +from openpmd_viewer.openpmd_timeseries import FieldMetaInformation + + +class OpenPMDDataReader(DataReader): + """ + Class perfroming the access to openPMD data files. + + This class is derived from the original openPMD `DataReader` to extend its + functionality and provide a method which returns only the field metadata. + """ + + def __init__(self, backend): + """ Initialize class. """ + super().__init__(backend) + + def read_field_metadata(self, iteration, field_name, component_name): + """ + Read the field metadata. + + Parameters: + ----------- + iteration : int + The iteration at which the parameters should be extracted. + field_name : str + Name of the field (e.g., `'E'`, `'B'`, `'rho'`, etc.). + component_name : str + Name of the field component (e.g., `'r'`, `'x'`, `'t'`, etc.) + + Returns: + -------- + A dictionary with the time, field and axis metadata. + """ + # Initialize metadata dictionary. + md = {} + + # Read basic openPMD parameters. + t, params = self.read_openPMD_params(iteration) + + # Time metadata. + md['time'] = {} + md['time']['value'] = t + md['time']['units'] = 's' + + # Field metadata. + md['field'] = {} + field_geometry = params['fields_metadata'][field_name]['geometry'] + axis_labels = params['fields_metadata'][field_name]['axis_labels'] + field_units = determine_field_units(field_name) + md['field']['geometry'] = field_geometry + md['field']['axis_labels'] = axis_labels + md['field']['units'] = field_units + + # Axis metadata. + md['axis'] = {} + info = self.get_field_meta_info( + iteration, field_name, component_name, axis_labels, field_geometry) + for axis in axis_labels: + md['axis'][axis] = {} + md['axis'][axis]['units'] = 'm' + md['axis'][axis]['array'] = getattr(info, axis) + md['axis'][axis]['spacing'] = getattr(info, 'd'+axis) + md['axis'][axis]['min'] = getattr(info, axis+'min') + md['axis'][axis]['max'] = getattr(info, axis+'max') + + return md + + def get_field_meta_info(self, iteration, field, comp, axis_labels, + geometry): + """ Get the `FieldMetaInformation` of the field. """ + if self.backend == 'h5py': + filename = self.iteration_to_file[iteration] + if geometry in ['thetaMode']: + return read_circ_field_metadata_h5py( + filename, iteration, field, comp) + elif geometry in ["1dcartesian", "2dcartesian", "3dcartesian"]: + return read_cartesian_field_metadata_h5py( + filename, iteration, field, comp, axis_labels) + elif self.backend == 'openpmd-api': + if geometry in ['thetaMode']: + return read_circ_field_metadata_io( + self.series, iteration, field, comp) + elif geometry in ["1dcartesian", "2dcartesian", "3dcartesian"]: + return read_cartesian_field_metadata_io( + self.series, iteration, field, comp) + + +def read_cartesian_field_metadata_h5py(filename, iteration, field_name, + component_name, axis_labels): + """ + Get `FieldMetaInformation` of cartesian field using `h5py` backend. + + Parameters: + ----------- + filename : str + The absolute path to the HDF5 file. + iteration : int + The iteration at which to obtain the data. + field_name : string + Which field to extract. + component_name : string, optional + Which component of the field to extract. + axis_labels: list of strings + The name of the dimensions of the array (e.g. ['x', 'y', 'z']). + """ + # Open the HDF5 file + dfile = h5py.File(filename, 'r') + # Extract the dataset and and corresponding group + if component_name is None: + field_path = field_name + else: + field_path = fr.join_infile_path(field_name, component_name) + group, dset = fr.find_dataset(dfile, iteration, field_path) + + # Extract the metainformation + shape = list(fr.get_shape(dset)) + axes = {i: axis_labels[i] for i in range(len(axis_labels))} + info = FieldMetaInformation( + axes, shape, group.attrs['gridSpacing'], + group.attrs['gridGlobalOffset'], group.attrs['gridUnitSI'], + dset.attrs['position']) + return info + + +def read_cartesian_field_metadata_io(series, iteration, field_name, + component_name, axis_labels): + """ + Get `FieldMetaInformation` of cartesian field using `io` backend. + + Parameters: + ----------- + series : openpmd_api.Series + An open, readable openPMD-api series object. + iteration : int + The iteration at which to obtain the data. + field_name : string + Which field to extract. + component_name : string, optional + Which component of the field to extract. + axis_labels: list of strings + The name of the dimensions of the array (e.g. ['x', 'y', 'z']). + """ + it = series.iterations[iteration] + + # Extract the dataset and and corresponding group + field = it.meshes[field_name] + if field.scalar: + component = next(field.items())[1] + else: + component = field[component_name] + + # Extract the metainformation + shape = component.shape + axes = {i: axis_labels[i] for i in range(len(axis_labels))} + info = FieldMetaInformation( + axes, shape, + field.grid_spacing, field.grid_global_offset, + field.grid_unit_SI, component.position) + return info + + +def read_circ_field_metadata_h5py(filename, iteration, field_name, + component_name): + """ + Get `FieldMetaInformation` of thetaMode field using `h5py` backend. + + Parameters: + ----------- + filename : str + The absolute path to the HDF5 file. + iteration : int + The iteration at which to obtain the data. + field_name : string + Which field to extract. + component_name : string, optional + Which component of the field to extract. + """ + # Open the HDF5 file + dfile = h5py.File(filename, 'r') + # Extract the dataset and and corresponding group + if component_name is None: + field_path = field_name + else: + field_path = fr.join_infile_path(field_name, component_name) + group, dset = fr.find_dataset(dfile, iteration, field_path) + + # Extract the metainformation + Nm, Nr, Nz = fr.get_shape(dset) + info = FieldMetaInformation( + {0: 'r', 1: 'z'}, (Nr, Nz), + group.attrs['gridSpacing'], group.attrs['gridGlobalOffset'], + group.attrs['gridUnitSI'], dset.attrs['position'], thetaMode=True) + + return info + + +def read_circ_field_metadata_io(series, iteration, field_name, component_name): + """ + Get `FieldMetaInformation` of thetaMode field using `io` backend. + + Parameters: + ----------- + series : openpmd_api.Series + An open, readable openPMD-api series object. + iteration : int + The iteration at which to obtain the data. + field_name : string + Which field to extract. + component_name : string, optional + Which component of the field to extract. + """ + it = series.iterations[iteration] + + # Extract the dataset and and corresponding group + field = it.meshes[field_name] + if field.scalar: + component = next(field.items())[1] + else: + component = field[component_name] + + # Extract the metainformation + Nm, Nr, Nz = component.shape + info = FieldMetaInformation( + {0: 'r', 1: 'z'}, (Nr, Nz), + field.grid_spacing, field.grid_global_offset, + field.grid_unit_SI, component.position, thetaMode=True) + return info + + +def determine_field_units(field_name): + """ Return the corresponding units of the field. """ + # TODO: Make more robust implementation using unit_dimension attributes. + if field_name == 'E': + return 'V/m' + elif field_name == 'B': + return 'T' + elif field_name == 'rho': + return 'C/m^3' + elif field_name == 'J': + return 'A' + else: + return '' From 54643b92c57ee0b1b2f1d526fbcbbc938aaf7ced Mon Sep 17 00:00:00 2001 From: AngelFP Date: Wed, 18 Aug 2021 13:44:38 +0200 Subject: [PATCH 2/4] Use new `OpenPMDDataReader` class --- visualpic/data_reading/field_readers.py | 43 +---------------------- visualpic/data_reading/folder_scanners.py | 4 +-- 2 files changed, 3 insertions(+), 44 deletions(-) diff --git a/visualpic/data_reading/field_readers.py b/visualpic/data_reading/field_readers.py index 7304871..303c081 100644 --- a/visualpic/data_reading/field_readers.py +++ b/visualpic/data_reading/field_readers.py @@ -479,45 +479,4 @@ def _read_field_metadata(self, file_path, iteration, field_path): field, *comp = field_path.split('/') if len(comp) > 0: comp = comp[0] - md = {} - t, params = self._opmd_reader.read_openPMD_params(iteration) - md['time'] = {} - md['time']['value'] = t - md['time']['units'] = 's' - field_geometry = params['fields_metadata'][field]['geometry'] - axis_labels = params['fields_metadata'][field]['axis_labels'] - field_units = self._determine_field_units(field) - md['field'] = {} - md['field']['units'] = field_units - md['field']['geometry'] = field_geometry - md['field']['axis_labels'] = axis_labels - ax_el, ax_lims = self._opmd_reader.get_grid_parameters( - iteration, [field], params['fields_metadata']) - axes = ax_el.keys() - md['axis'] = {} - for axis in axes: - md['axis'][axis] = {} - md['axis'][axis]['units'] = 'm' - ax_min = ax_lims[axis][0] - ax_max = ax_lims[axis][1] - ax_els = ax_el[axis] - if field_geometry in ['cylindrical', 'thetaMode'] and axis == 'r': - ax_min = -ax_max - ax_els += ax_el[axis] - # FIXME this does not differentiate between - # node-centered / cell-centered fields. FieldMetaInformation - # does it properly - md['axis'][axis]['array'] = np.linspace(ax_min, ax_max, ax_els) - return md - - def _determine_field_units(self, field): - if field == 'E': - return 'V/m' - elif field == 'B': - return 'T' - elif field == 'rho': - return 'C/m^3' - elif field == 'J': - return 'A' - else: - return '' + return self._opmd_reader.read_field_metadata(iteration, field, comp) diff --git a/visualpic/data_reading/folder_scanners.py b/visualpic/data_reading/folder_scanners.py index b19cef5..55a8e2c 100644 --- a/visualpic/data_reading/folder_scanners.py +++ b/visualpic/data_reading/folder_scanners.py @@ -12,7 +12,7 @@ import numpy as np from h5py import File as H5F -from openpmd_viewer.openpmd_timeseries.data_reader import DataReader +from .openpmd_data_reader import OpenPMDDataReader import visualpic.data_reading.field_readers as fr import visualpic.data_reading.particle_readers as pr @@ -77,7 +77,7 @@ def __init__(self, opmd_backend='h5py'): Possible values are 'h5py' or 'openpmd-api'. """ - self.opmd_reader = DataReader(opmd_backend) + self.opmd_reader = OpenPMDDataReader(opmd_backend) self.field_reader = fr.OpenPMDFieldReader(self.opmd_reader) self.particle_reader = pr.OpenPMDParticleReader(self.opmd_reader) self.unit_converter = uc.OpenPMDUnitConverter() From 51c36f58d9144874b2fc1a7d972bf4f59b59606e Mon Sep 17 00:00:00 2001 From: AngelFP Date: Wed, 18 Aug 2021 14:00:45 +0200 Subject: [PATCH 3/4] Fix bug when field is scalar --- visualpic/data_reading/field_readers.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/visualpic/data_reading/field_readers.py b/visualpic/data_reading/field_readers.py index 303c081..dfb203f 100644 --- a/visualpic/data_reading/field_readers.py +++ b/visualpic/data_reading/field_readers.py @@ -55,19 +55,24 @@ def _readjust_metadata(self, field_metadata, slice_dir_i, slice_dir_j, geom = field_metadata['field']['geometry'] if geom in ['cylindrical', 'thetaMode'] and theta is None: r_md = field_metadata['axis']['r'] + z_md = field_metadata['axis']['z'] # Check if resolution should be reduced if max_resolution_3d is not None: - z = field_metadata['axis']['z']['array'] + z = z_md['array'] r = r_md['array'] nr = len(r) - 1 nz = len(z) - 1 max_res_lon, max_res_transv = max_resolution_3d if nz > max_res_lon: excess_z = int(np.round(nz/max_res_lon)) - field_metadata['axis']['z']['array'] = z[::excess_z] + z_md['array'] = z[::excess_z] + z_md['min'] = z_md['array'][0] + z_md['max'] = z_md['array'][-1] if nr > max_res_transv: excess_r = int(np.round(nr/max_res_transv)) r_md['array'] = r[::excess_r] + r_md['min'] = z_md['array'][0] + r_md['max'] = z_md['array'][-1] # if nr > max_res_transv: # r = zoom(r, max_res_transv/nr, order=1) # r_md['array'] = r @@ -476,7 +481,11 @@ def _read_field_theta( return fld def _read_field_metadata(self, file_path, iteration, field_path): + # Get name of field and component. field, *comp = field_path.split('/') if len(comp) > 0: comp = comp[0] + else: + comp = None + # Read metadata. return self._opmd_reader.read_field_metadata(iteration, field, comp) From 00b2849da76d16113ec17509a9897ddad7e461f4 Mon Sep 17 00:00:00 2001 From: AngelFP Date: Wed, 18 Aug 2021 14:05:30 +0200 Subject: [PATCH 4/4] Add module header --- visualpic/data_reading/openpmd_data_reader.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/visualpic/data_reading/openpmd_data_reader.py b/visualpic/data_reading/openpmd_data_reader.py index 9e47d8f..871bda4 100644 --- a/visualpic/data_reading/openpmd_data_reader.py +++ b/visualpic/data_reading/openpmd_data_reader.py @@ -1,3 +1,12 @@ +""" +This module contains the `OpenPMDDataReader`, a class derived from the +openPMD-viewer `DataReader` which has been expanded with a method for returning +only the field metadata. + +Some of the methods below include code from openPMD-viewer (see at +https://github.com/openPMD/openPMD-viewer). +""" + import h5py from openpmd_viewer.openpmd_timeseries.data_reader import DataReader from openpmd_viewer.openpmd_timeseries.data_reader.h5py_reader import (