From 6603ae4d675e58327856047cc9f804e140ad8724 Mon Sep 17 00:00:00 2001 From: Alberto Martinez de la Ossa Date: Wed, 10 Jun 2020 17:28:50 +0200 Subject: [PATCH 01/16] updated comment --- visualpic/data_reading/particle_readers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/visualpic/data_reading/particle_readers.py b/visualpic/data_reading/particle_readers.py index 56027ea..2b0edd1 100644 --- a/visualpic/data_reading/particle_readers.py +++ b/visualpic/data_reading/particle_readers.py @@ -116,7 +116,7 @@ def _read_component_data(self, file_path, species, component): b = data[:, 1] data = 1/2*(a+b)*(a+b+1)+b """ - # This block is inconsistent with how is done in OSIRIS + # This block is inconsistent with how is done in OSIRIS (above) if component == 'q': n_cells = file_handle.attrs['NX'] sim_size = (file_handle.attrs['XMAX'] - file_handle.attrs['XMIN']) From 15b03f7528ee8b949a2ffc58bdf4996e4a6dc541 Mon Sep 17 00:00:00 2001 From: Alberto Martinez de la Ossa Date: Wed, 10 Jun 2020 18:25:09 +0200 Subject: [PATCH 02/16] reapplying changes by Alberto --- .../derived_field_definitions.py | 9 ++-- visualpic/data_handling/unit_converters.py | 45 ++++++++++++++++--- visualpic/data_reading/field_readers.py | 4 +- visualpic/data_reading/particle_readers.py | 25 ++++++----- 4 files changed, 61 insertions(+), 22 deletions(-) diff --git a/visualpic/data_handling/derived_field_definitions.py b/visualpic/data_handling/derived_field_definitions.py index 15614d8..4c657d4 100644 --- a/visualpic/data_handling/derived_field_definitions.py +++ b/visualpic/data_handling/derived_field_definitions.py @@ -125,11 +125,14 @@ def calculate_norm_vector_pot(data_list, sim_geometry, sim_params): elif sim_geometry == 'thetaMode': Ez, Er, Et = data_list E2 = Ez**2 + Er**2 + Et**2 - return np.sqrt(7.3e-11 * l_0**2 * ct.c * ct.epsilon_0 / 2 * E2) - + # return np.sqrt(7.3e-11 * l_0**2 * ct.c * ct.epsilon_0 / 2 * E2) + k_0 = 2. * np.pi / l_0 + # return (ct.e / ct.m_e / ct.c**2) * np.sqrt(E2) / k_0 + return np.sqrt(E2) / k_0 norm_vector_pot = {'name': 'a', - 'units': 'm_e*c^2/e', + # 'units': 'm_e*c^2/e', + 'units': 'V', 'requirements': {'1d': ['Ez'], '2dcartesian': ['Ez', 'Ex'], '3dcartesian': ['Ez', 'Ex', 'Ey'], diff --git a/visualpic/data_handling/unit_converters.py b/visualpic/data_handling/unit_converters.py index 6c53791..505ed60 100644 --- a/visualpic/data_handling/unit_converters.py +++ b/visualpic/data_handling/unit_converters.py @@ -13,7 +13,8 @@ import aptools.plasma_accel.general_equations as ge -length_conversion = {'km': 1e-3, +length_conversion = {'m': 1., + 'km': 1e-3, 'mm': 1e3, 'um': 1e6, 'nm': 1e9} @@ -64,6 +65,7 @@ intensity_conversion = {'W/cm^2': 1e-4} +potential_conversion = {'V': 1} class UnitConverter(): def __init__(self): @@ -76,7 +78,8 @@ def __init__(self): 'V/m^2': efieldgradient_conversion, 'T/m': bfieldgradient_conversion, 'm_e*c': momentum_conversion, - 'W/m^2': intensity_conversion} + 'W/m^2': intensity_conversion, + 'V': potential_conversion} self.si_units = list(self.conversion_factors.keys()) @@ -243,8 +246,11 @@ def __init__(self, plasma_density=None): '1/\\omega_p': [1./w_p, 's'], 'c/\\omega_p': [ct.c/w_p, 'm'], 'm_ec\\omega_pe^{-1}': [E_0, 'V/m'], - 'e\\omega_p^3/c^3': [(w_p/ct.c)**3, '1/m^3'], - 'm_ec': [1., 'm_e*c']} + # 'e\\omega_p^3/c^3': [(w_p/ct.c)**3, '1/m^3'], + 'e\\omega_p^3/c^3': [ct.e * self.plasma_density, 'C/m^3'], + 'm_ec': [1., 'm_e*c'], + 'm_e*c^2/e': [ct.m_e*ct.c**2/ct.e, 'V'] + } else: self.osiris_unit_conversion = None super().__init__() @@ -272,7 +278,34 @@ def convert_data_to_si(self, data, data_units, metadata=None): class HiPACEUnitConverter(UnitConverter): def __init__(self, plasma_density=None): self.plasma_density = plasma_density + if plasma_density is not None: + w_p = ge.plasma_frequency(plasma_density*1e-6) + E_0 = ge.plasma_cold_non_relativisct_wave_breaking_field( + plasma_density*1e-6) + self.hipace_unit_conversion = {'1/\\omega_p': [1/w_p, 's'], + 'c/\\omega_p': [ct.c/w_p, 'm'], + 'E_0': [E_0, 'V/m'], + 'n_0': [ct.e * self.plasma_density, 'C/m^3'], + 'm_ec': [ct.m_e*ct.c, 'J*s/m']} + else: + self.hipace_unit_conversion = None + super().__init__() def convert_data_to_si(self, data, data_units, metadata=None): - raise NotImplementedError( - 'HiPACE unit conversion not yet implemented.') + if self.hipace_unit_conversion is not None: + if data_units in self.hipace_unit_conversion: + conv_factor, si_units = self.hipace_unit_conversion[data_units] + elif data_units == 'qnorm': + n_cells = metadata['grid']['resolution'] + sim_size = metadata['grid']['size'] + cell_vol = np.prod(sim_size/n_cells) + s_d = ge.plasma_skin_depth(self.plasma_density*1e-6) + conv_factor = cell_vol * self.plasma_density * s_d**3 * ct.e + si_units = 'C' + else: + raise ValueError('Unsupported units: {}.'.format(data_units)) + return data*conv_factor, si_units + + else: + raise ValueError('Could not perform unit conversion.' + ' Plasma density value not provided.') diff --git a/visualpic/data_reading/field_readers.py b/visualpic/data_reading/field_readers.py index fa5eb5d..1188d31 100644 --- a/visualpic/data_reading/field_readers.py +++ b/visualpic/data_reading/field_readers.py @@ -274,9 +274,9 @@ def _get_field_units(self, file_path): # HiPACE does not store unit information in data files file = os.path.split(file_path)[-1] if 'field' in file: - return 'm_e c \\omega_p e^{-1}' + return 'E_0' elif 'density' in file: - 'e \\omega_p^3/ c^3' + return 'n_0' else: return '' diff --git a/visualpic/data_reading/particle_readers.py b/visualpic/data_reading/particle_readers.py index 58db0b0..2b0edd1 100644 --- a/visualpic/data_reading/particle_readers.py +++ b/visualpic/data_reading/particle_readers.py @@ -75,9 +75,12 @@ def _read_component_metadata(self, file_path, species, component): metadata['time']['units'] = self._numpy_bytes_to_string( file_handle.attrs['TIME UNITS'][0]) metadata['grid'] = {} - metadata['grid']['resolution'] = file_handle.attrs['NX'] - max_range = file_handle.attrs['XMAX'] - min_range = file_handle.attrs['XMIN'] + simdata_path = '/SIMULATION' + if simdata_path not in file_handle.keys(): + simdata_path = '/' + metadata['grid']['resolution'] = file_handle[simdata_path].attrs['NX'] + max_range = file_handle[simdata_path].attrs['XMAX'] + min_range = file_handle[simdata_path].attrs['XMIN'] metadata['grid']['size'] = max_range - min_range grid_range = [] for x_min, x_max in zip(min_range, max_range): @@ -113,28 +116,28 @@ def _read_component_data(self, file_path, species, component): b = data[:, 1] data = 1/2*(a+b)*(a+b+1)+b """ - # This block is inconsistent with how is done in OsirisParticleReader + # This block is inconsistent with how is done in OSIRIS (above) if component == 'q': n_cells = file_handle.attrs['NX'] sim_size = (file_handle.attrs['XMAX'] - file_handle.attrs['XMIN']) cell_vol = np.prod(sim_size/n_cells) data *= cell_vol - return np.array(data) """ - + return np.array(data) + def _read_component_metadata(self, file_path, species, component): metadata = {} with H5F(file_path, 'r') as file_handle: if component in ['x', 'y', 'z']: - units = 'c/ \omega_p' + units = 'c/\\omega_p' elif component in ['px', 'py', 'pz']: - units = 'm_e c' + units = 'm_ec' elif component == 'q': - units = 'e n_p c^3 / \\omega_p^3' + units = 'qnorm' metadata['units'] = units metadata['time'] = {} metadata['time']['value'] = file_handle.attrs['TIME'][0] - metadata['time']['units'] = '1/ \\omega_p' + metadata['time']['units'] = '1/\\omega_p' metadata['grid'] = {} metadata['grid']['resolution'] = file_handle.attrs['NX'] max_range = file_handle.attrs['XMAX'] @@ -144,7 +147,7 @@ def _read_component_metadata(self, file_path, species, component): for x_min, x_max in zip(min_range, max_range): grid_range.append([x_min, x_max]) metadata['grid']['range'] = grid_range - metadata['grid']['size_units'] = '\\omega_p/c' + metadata['grid']['size_units'] = 'c/\\omega_p' return metadata From 48b39e24f8cccf591cfe1ad7ea1625c08e368fc4 Mon Sep 17 00:00:00 2001 From: AngelFP Date: Mon, 15 Jun 2020 13:32:30 +0200 Subject: [PATCH 03/16] Change potential conversion --- visualpic/data_handling/unit_converters.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/visualpic/data_handling/unit_converters.py b/visualpic/data_handling/unit_converters.py index 505ed60..41c50d5 100644 --- a/visualpic/data_handling/unit_converters.py +++ b/visualpic/data_handling/unit_converters.py @@ -65,7 +65,10 @@ intensity_conversion = {'W/cm^2': 1e-4} -potential_conversion = {'V': 1} + +potential_conversion = {'kV': 1e-3, + 'MV': 1e-6} + class UnitConverter(): def __init__(self): From 2440cdccc4dc31d10809426d70cdecf624f54dc0 Mon Sep 17 00:00:00 2001 From: AngelFP Date: Mon, 15 Jun 2020 13:33:32 +0200 Subject: [PATCH 04/16] Remove unnecessary conversion --- visualpic/data_handling/unit_converters.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/visualpic/data_handling/unit_converters.py b/visualpic/data_handling/unit_converters.py index 41c50d5..01bab1e 100644 --- a/visualpic/data_handling/unit_converters.py +++ b/visualpic/data_handling/unit_converters.py @@ -13,8 +13,7 @@ import aptools.plasma_accel.general_equations as ge -length_conversion = {'m': 1., - 'km': 1e-3, +length_conversion = {'km': 1e-3, 'mm': 1e3, 'um': 1e6, 'nm': 1e9} From 8aa05cff2d15a9db5a3e3b937a51e0687ffb0d60 Mon Sep 17 00:00:00 2001 From: AngelFP Date: Mon, 15 Jun 2020 13:52:13 +0200 Subject: [PATCH 05/16] Fix long line --- visualpic/data_handling/unit_converters.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/visualpic/data_handling/unit_converters.py b/visualpic/data_handling/unit_converters.py index 01bab1e..add949a 100644 --- a/visualpic/data_handling/unit_converters.py +++ b/visualpic/data_handling/unit_converters.py @@ -284,11 +284,13 @@ def __init__(self, plasma_density=None): w_p = ge.plasma_frequency(plasma_density*1e-6) E_0 = ge.plasma_cold_non_relativisct_wave_breaking_field( plasma_density*1e-6) - self.hipace_unit_conversion = {'1/\\omega_p': [1/w_p, 's'], - 'c/\\omega_p': [ct.c/w_p, 'm'], - 'E_0': [E_0, 'V/m'], - 'n_0': [ct.e * self.plasma_density, 'C/m^3'], - 'm_ec': [ct.m_e*ct.c, 'J*s/m']} + self.hipace_unit_conversion = { + '1/\\omega_p': [1/w_p, 's'], + 'c/\\omega_p': [ct.c/w_p, 'm'], + 'E_0': [E_0, 'V/m'], + 'n_0': [ct.e * self.plasma_density, 'C/m^3'], + 'm_ec': [ct.m_e*ct.c, 'J*s/m'] + } else: self.hipace_unit_conversion = None super().__init__() From 204f96d83e6fef0aa6ae0022196a9610af2311e4 Mon Sep 17 00:00:00 2001 From: AngelFP Date: Mon, 15 Jun 2020 13:57:57 +0200 Subject: [PATCH 06/16] Add missing 'axis_labels' metadata --- visualpic/data_reading/field_readers.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/visualpic/data_reading/field_readers.py b/visualpic/data_reading/field_readers.py index 1188d31..91b871b 100644 --- a/visualpic/data_reading/field_readers.py +++ b/visualpic/data_reading/field_readers.py @@ -162,6 +162,7 @@ def _read_field_metadata(self, file_path, field_path): md['field'] = {} md['field']['units'] = field_units md['field']['geometry'] = field_geometry + # TODO: check correct order of labels if field_geometry == "3dcartesian": axis_labels = ['x', 'y', 'z'] elif field_geometry == "2dcartesian": @@ -264,6 +265,8 @@ def _read_field_metadata(self, file_path, field_path): md['field'] = {} md['field']['units'] = field_units md['field']['geometry'] = '3dcartesian' + # TODO: check correct order of labels + md['field']['axis_labels'] = ['x', 'y', 'z'] md['axis'] = self._get_axis_data(file, field_shape) md['time'] = self._get_time_data(file) file.close() From f089368264845cc451ca731d39524daa2e6fe0b2 Mon Sep 17 00:00:00 2001 From: Alberto Martinez de la Ossa Date: Mon, 15 Jun 2020 17:48:16 +0200 Subject: [PATCH 07/16] fix an inconsistency on the axis units of HiPACE --- visualpic/data_reading/field_readers.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/visualpic/data_reading/field_readers.py b/visualpic/data_reading/field_readers.py index 91b871b..64bf53b 100644 --- a/visualpic/data_reading/field_readers.py +++ b/visualpic/data_reading/field_readers.py @@ -291,17 +291,17 @@ def _get_axis_data(self, file, field_shape): """ Returns dictionary with the array and units of each field axis """ axis_data = {} axis_data['z'] = {} - axis_data["z"]["units"] = 'c/ \\omega_p' + axis_data["z"]["units"] = 'c/\\omega_p' axis_data["z"]["array"] = np.linspace(file.attrs['XMIN'][0], file.attrs['XMAX'][0], field_shape[0]+1) axis_data['x'] = {} - axis_data["x"]["units"] = 'c/ \\omega_p' + axis_data["x"]["units"] = 'c/\\omega_p' axis_data["x"]["array"] = np.linspace(file.attrs['XMIN'][1], file.attrs['XMAX'][1], field_shape[1]+1) axis_data['y'] = {} - axis_data["y"]["units"] = 'c/ \\omega_p' + axis_data["y"]["units"] = 'c/\\omega_p' axis_data["y"]["array"] = np.linspace(file.attrs['XMIN'][2], file.attrs['XMAX'][2], field_shape[2]+1) @@ -311,7 +311,7 @@ def _get_time_data(self, file): """ Returns dictionary with value and units of the simulation time """ time_data = {} time_data["value"] = file.attrs["TIME"][0] - time_data["units"] = '1/ \\omega_p' + time_data["units"] = '1/\\omega_p' return time_data From ff4e1a9e71e9ceece3e148d8e7c17197b792e10e Mon Sep 17 00:00:00 2001 From: AngelFP Date: Tue, 16 Jun 2020 09:57:22 +0200 Subject: [PATCH 08/16] Implement separate 'a' and 'A' fields --- .../derived_field_definitions.py | 29 +++++++++++++++---- 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/visualpic/data_handling/derived_field_definitions.py b/visualpic/data_handling/derived_field_definitions.py index 4c657d4..fec3741 100644 --- a/visualpic/data_handling/derived_field_definitions.py +++ b/visualpic/data_handling/derived_field_definitions.py @@ -108,8 +108,8 @@ def calculate_intensity(data_list, sim_geometry, sim_params): derived_field_definitions.append(intensity) -# Normalized vector potential -def calculate_norm_vector_pot(data_list, sim_geometry, sim_params): +# Vector potential +def calculate_vector_pot(data_list, sim_geometry, sim_params): l_0 = sim_params['lambda_0'] if sim_geometry == '1d': Ez = data_list[0] @@ -125,14 +125,31 @@ def calculate_norm_vector_pot(data_list, sim_geometry, sim_params): elif sim_geometry == 'thetaMode': Ez, Er, Et = data_list E2 = Ez**2 + Er**2 + Et**2 - # return np.sqrt(7.3e-11 * l_0**2 * ct.c * ct.epsilon_0 / 2 * E2) k_0 = 2. * np.pi / l_0 - # return (ct.e / ct.m_e / ct.c**2) * np.sqrt(E2) / k_0 return np.sqrt(E2) / k_0 + +vector_pot = {'name': 'A', + 'units': 'V', + 'requirements': {'1d': ['Ez'], + '2dcartesian': ['Ez', 'Ex'], + '3dcartesian': ['Ez', 'Ex', 'Ey'], + '2dcylindrical': [], + 'thetaMode': ['Ez', 'Er', 'Et']}, + 'recipe': calculate_vector_pot} + + +derived_field_definitions.append(vector_pot) + + +# Normalized vector potential +def calculate_norm_vector_pot(data_list, sim_geometry, sim_params): + A = calculate_vector_pot(data_list, sim_geometry, sim_params) + return A / ((ct.m_e * ct.c**2) / ct.e) + + norm_vector_pot = {'name': 'a', - # 'units': 'm_e*c^2/e', - 'units': 'V', + 'units': 'm_e*c^2/e', 'requirements': {'1d': ['Ez'], '2dcartesian': ['Ez', 'Ex'], '3dcartesian': ['Ez', 'Ex', 'Ey'], From 0d1c860cfdc1f5593a191f68f71e2205fb42fed6 Mon Sep 17 00:00:00 2001 From: AngelFP Date: Tue, 16 Jun 2020 09:57:54 +0200 Subject: [PATCH 09/16] Add unit conversion for vector potentials --- visualpic/data_handling/unit_converters.py | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/visualpic/data_handling/unit_converters.py b/visualpic/data_handling/unit_converters.py index add949a..29e4678 100644 --- a/visualpic/data_handling/unit_converters.py +++ b/visualpic/data_handling/unit_converters.py @@ -65,12 +65,22 @@ intensity_conversion = {'W/cm^2': 1e-4} -potential_conversion = {'kV': 1e-3, +potential_conversion = {'m_e*c^2/e': 1 / ((ct.m_e * ct.c**2) / ct.e), + 'kV': 1e-3, 'MV': 1e-6} +norm_potential_conversion = {'V': ct.m_e * ct.c**2 / ct.e, + 'kV': ct.m_e * ct.c**2 / ct.e * 1e-3, + 'MV': ct.m_e * ct.c**2 / ct.e * 1e-6} + + class UnitConverter(): def __init__(self): + """ Initialize unit converter. """ + # For convenience, due to their common use, the momentum units 'm_e*c' + # and the normalized potential units 'm_e*c^2/e' are also considered SI + # units. self.conversion_factors = {'m': length_conversion, 's': time_conversion, 'rad': angle_conversion, @@ -81,7 +91,8 @@ def __init__(self): 'T/m': bfieldgradient_conversion, 'm_e*c': momentum_conversion, 'W/m^2': intensity_conversion, - 'V': potential_conversion} + 'V': potential_conversion, + 'm_e*c^2/e': norm_potential_conversion} self.si_units = list(self.conversion_factors.keys()) @@ -113,8 +124,9 @@ def convert_field_units(self, field_data, field_md, axes_to_convert, convert_time) # convert field data to desired units + field_units = field_md['field']['units'] if convert_field and (target_field_units != 'SI' and - target_field_units not in self.si_units): + target_field_units != field_units): field_units = field_md['field']['units'] field_data = self.convert_data(field_data, field_units, target_field_units) @@ -250,8 +262,7 @@ def __init__(self, plasma_density=None): 'm_ec\\omega_pe^{-1}': [E_0, 'V/m'], # 'e\\omega_p^3/c^3': [(w_p/ct.c)**3, '1/m^3'], 'e\\omega_p^3/c^3': [ct.e * self.plasma_density, 'C/m^3'], - 'm_ec': [1., 'm_e*c'], - 'm_e*c^2/e': [ct.m_e*ct.c**2/ct.e, 'V'] + 'm_ec': [1., 'm_e*c'] } else: self.osiris_unit_conversion = None From 6d623888cf7c22fb03f3696ababbe0fc054a0b2e Mon Sep 17 00:00:00 2001 From: AngelFP Date: Tue, 16 Jun 2020 10:07:05 +0200 Subject: [PATCH 10/16] Remove commented block --- visualpic/data_reading/particle_readers.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/visualpic/data_reading/particle_readers.py b/visualpic/data_reading/particle_readers.py index 2b0edd1..6d02af4 100644 --- a/visualpic/data_reading/particle_readers.py +++ b/visualpic/data_reading/particle_readers.py @@ -115,14 +115,6 @@ def _read_component_data(self, file_path, species, component): a = data[:, 0] b = data[:, 1] data = 1/2*(a+b)*(a+b+1)+b - """ - # This block is inconsistent with how is done in OSIRIS (above) - if component == 'q': - n_cells = file_handle.attrs['NX'] - sim_size = (file_handle.attrs['XMAX'] - file_handle.attrs['XMIN']) - cell_vol = np.prod(sim_size/n_cells) - data *= cell_vol - """ return np.array(data) def _read_component_metadata(self, file_path, species, component): From 9e39ff3a0bd3eee85413671161ef68f3d3250ccd Mon Sep 17 00:00:00 2001 From: AngelFP Date: Tue, 16 Jun 2020 10:09:39 +0200 Subject: [PATCH 11/16] Remove duplicate line --- visualpic/data_handling/unit_converters.py | 1 - 1 file changed, 1 deletion(-) diff --git a/visualpic/data_handling/unit_converters.py b/visualpic/data_handling/unit_converters.py index 29e4678..3991b6b 100644 --- a/visualpic/data_handling/unit_converters.py +++ b/visualpic/data_handling/unit_converters.py @@ -127,7 +127,6 @@ def convert_field_units(self, field_data, field_md, field_units = field_md['field']['units'] if convert_field and (target_field_units != 'SI' and target_field_units != field_units): - field_units = field_md['field']['units'] field_data = self.convert_data(field_data, field_units, target_field_units) field_md['field']['units'] = target_field_units From 200cdf2bcb49ebbb6684bb618189e3b775b53bf2 Mon Sep 17 00:00:00 2001 From: Alberto Martinez de la Ossa Date: Tue, 16 Jun 2020 13:08:25 +0200 Subject: [PATCH 12/16] implemented fix for accessing simulation parameters in newer versions of Osiris --- visualpic/data_reading/field_readers.py | 20 +++++++++++++------- visualpic/data_reading/particle_readers.py | 9 +++++++-- 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/visualpic/data_reading/field_readers.py b/visualpic/data_reading/field_readers.py index 64bf53b..17adf8e 100644 --- a/visualpic/data_reading/field_readers.py +++ b/visualpic/data_reading/field_readers.py @@ -181,7 +181,10 @@ def _read_field_metadata(self, file_path, field_path): def _get_field_units(self, file, field_path): """ Returns the field units""" - return self._numpy_bytes_to_string(file[field_path].attrs["UNITS"][0]) + if "UNITS" in file[field_path].attrs: + return self._numpy_bytes_to_string(file[field_path].attrs["UNITS"][0]) + else: + return self._numpy_bytes_to_string(file.attrs["UNITS"][0]) def _get_field_shape(self, file, field_path): """ Returns shape of field array""" @@ -200,24 +203,27 @@ def _get_axis_data(self, file, field_path, field_geometry, field_shape): """ Returns dictionary with the array and units of each field axis """ axis_data = {} axis_data['z'] = {} + simdata_path = '/SIMULATION' + if simdata_path not in file.keys(): + simdata_path = '/' axis_data["z"]["units"] = self._numpy_bytes_to_string( file['/AXIS/AXIS1'].attrs["UNITS"][0]) - axis_data["z"]["array"] = np.linspace(file.attrs['XMIN'][0], - file.attrs['XMAX'][0], + axis_data["z"]["array"] = np.linspace(file[simdata_path].attrs['XMIN'][0], + file[simdata_path].attrs['XMAX'][0], field_shape[-1]+1) if field_geometry in ["2dcartesian", "3dcartesian"]: axis_data['x'] = {} axis_data["x"]["units"] = self._numpy_bytes_to_string( file['/AXIS/AXIS2'].attrs["UNITS"][0]) - axis_data["x"]["array"] = np.linspace(file.attrs['XMIN'][1], - file.attrs['XMAX'][1], + axis_data["x"]["array"] = np.linspace(file[simdata_path].attrs['XMIN'][1], + file[simdata_path].attrs['XMAX'][1], field_shape[0]+1) if field_geometry == "3dcartesian": axis_data['y'] = {} axis_data["y"]["units"] = self._numpy_bytes_to_string( file['/AXIS/AXIS3'].attrs["UNITS"][0]) - axis_data["y"]["array"] = np.linspace(file.attrs['XMIN'][2], - file.attrs['XMAX'][2], + axis_data["y"]["array"] = np.linspace(file[simdata_path].attrs['XMIN'][2], + file[simdata_path].attrs['XMAX'][2], field_shape[1]+1) return axis_data diff --git a/visualpic/data_reading/particle_readers.py b/visualpic/data_reading/particle_readers.py index 2b0edd1..896a83c 100644 --- a/visualpic/data_reading/particle_readers.py +++ b/visualpic/data_reading/particle_readers.py @@ -68,8 +68,13 @@ def _read_component_metadata(self, file_path, species, component): metadata = {} with H5F(file_path, 'r') as file_handle: if component != 'tag': - metadata['units'] = self._numpy_bytes_to_string(file_handle[ - self.name_relations[component]].attrs['UNITS'][0]) + if 'QUANTS' in file_handle.attrs: + quantlist = [self._numpy_bytes_to_string(q) for q in file_handle.attrs['QUANTS']] + idx = quantlist.index(self.name_relations[component]) + metadata['units'] = self._numpy_bytes_to_string(file_handle.attrs['UNITS'][idx]) + else: + metadata['units'] = self._numpy_bytes_to_string(file_handle[self.name_relations[component]].attrs['UNITS'][0]) + metadata['time'] = {} metadata['time']['value'] = file_handle.attrs['TIME'][0] metadata['time']['units'] = self._numpy_bytes_to_string( From 60bad79f3584bdbdd8f21d310c342f5cbb5150f6 Mon Sep 17 00:00:00 2001 From: AngelFP Date: Tue, 16 Jun 2020 13:58:50 +0200 Subject: [PATCH 13/16] Fix formatting, add comments. --- visualpic/data_reading/field_readers.py | 25 ++++++++++--------- visualpic/data_reading/particle_readers.py | 28 +++++++++++++++------- 2 files changed, 33 insertions(+), 20 deletions(-) diff --git a/visualpic/data_reading/field_readers.py b/visualpic/data_reading/field_readers.py index 17adf8e..ed26127 100644 --- a/visualpic/data_reading/field_readers.py +++ b/visualpic/data_reading/field_readers.py @@ -181,10 +181,11 @@ def _read_field_metadata(self, file_path, field_path): def _get_field_units(self, file, field_path): """ Returns the field units""" + attr_path = '/' + # In older Osiris versions the units are in field_path. if "UNITS" in file[field_path].attrs: - return self._numpy_bytes_to_string(file[field_path].attrs["UNITS"][0]) - else: - return self._numpy_bytes_to_string(file.attrs["UNITS"][0]) + attr_path = field_path + return self._numpy_bytes_to_string(file[attr_path].attrs["UNITS"][0]) def _get_field_shape(self, file, field_path): """ Returns shape of field array""" @@ -201,29 +202,31 @@ def _determine_geometry(self, file): def _get_axis_data(self, file, field_path, field_geometry, field_shape): """ Returns dictionary with the array and units of each field axis """ - axis_data = {} - axis_data['z'] = {} simdata_path = '/SIMULATION' + # In older Osiris versions the simulation parameters are in '/'. if simdata_path not in file.keys(): simdata_path = '/' + sim_data = file[simdata_path] + axis_data = {} + axis_data['z'] = {} axis_data["z"]["units"] = self._numpy_bytes_to_string( file['/AXIS/AXIS1'].attrs["UNITS"][0]) - axis_data["z"]["array"] = np.linspace(file[simdata_path].attrs['XMIN'][0], - file[simdata_path].attrs['XMAX'][0], + axis_data["z"]["array"] = np.linspace(sim_data.attrs['XMIN'][0], + sim_data.attrs['XMAX'][0], field_shape[-1]+1) if field_geometry in ["2dcartesian", "3dcartesian"]: axis_data['x'] = {} axis_data["x"]["units"] = self._numpy_bytes_to_string( file['/AXIS/AXIS2'].attrs["UNITS"][0]) - axis_data["x"]["array"] = np.linspace(file[simdata_path].attrs['XMIN'][1], - file[simdata_path].attrs['XMAX'][1], + axis_data["x"]["array"] = np.linspace(sim_data.attrs['XMIN'][1], + sim_data.attrs['XMAX'][1], field_shape[0]+1) if field_geometry == "3dcartesian": axis_data['y'] = {} axis_data["y"]["units"] = self._numpy_bytes_to_string( file['/AXIS/AXIS3'].attrs["UNITS"][0]) - axis_data["y"]["array"] = np.linspace(file[simdata_path].attrs['XMIN'][2], - file[simdata_path].attrs['XMAX'][2], + axis_data["y"]["array"] = np.linspace(sim_data.attrs['XMIN'][2], + sim_data.attrs['XMAX'][2], field_shape[1]+1) return axis_data diff --git a/visualpic/data_reading/particle_readers.py b/visualpic/data_reading/particle_readers.py index 6e3eb1d..f2dfcda 100644 --- a/visualpic/data_reading/particle_readers.py +++ b/visualpic/data_reading/particle_readers.py @@ -67,25 +67,35 @@ def _read_component_data(self, file_path, species, component): def _read_component_metadata(self, file_path, species, component): metadata = {} with H5F(file_path, 'r') as file_handle: + # Read units. if component != 'tag': + osiris_name = self.name_relations[component] + # In new Osiris versions, the units are in a list in '/'. if 'QUANTS' in file_handle.attrs: - quantlist = [self._numpy_bytes_to_string(q) for q in file_handle.attrs['QUANTS']] - idx = quantlist.index(self.name_relations[component]) - metadata['units'] = self._numpy_bytes_to_string(file_handle.attrs['UNITS'][idx]) + units_path = '/' + quantlist = [self._numpy_bytes_to_string(q) + for q in file_handle.attrs['QUANTS']] + idx = quantlist.index(osiris_name) else: - metadata['units'] = self._numpy_bytes_to_string(file_handle[self.name_relations[component]].attrs['UNITS'][0]) - + units_path = osiris_name + idx = 0 + metadata['units'] = self._numpy_bytes_to_string( + file_handle[units_path].attrs['UNITS'][idx]) + # Read time data. metadata['time'] = {} metadata['time']['value'] = file_handle.attrs['TIME'][0] metadata['time']['units'] = self._numpy_bytes_to_string( file_handle.attrs['TIME UNITS'][0]) - metadata['grid'] = {} + # Read grid parameters. simdata_path = '/SIMULATION' + # In older Osiris versions the simulation parameters are in '/'. if simdata_path not in file_handle.keys(): simdata_path = '/' - metadata['grid']['resolution'] = file_handle[simdata_path].attrs['NX'] - max_range = file_handle[simdata_path].attrs['XMAX'] - min_range = file_handle[simdata_path].attrs['XMIN'] + sim_data = file_handle[simdata_path] + metadata['grid'] = {} + metadata['grid']['resolution'] = sim_data.attrs['NX'] + max_range = sim_data.attrs['XMAX'] + min_range = sim_data.attrs['XMIN'] metadata['grid']['size'] = max_range - min_range grid_range = [] for x_min, x_max in zip(min_range, max_range): From c2f16975ecf776e389b9aba24dd748e5c0c699da Mon Sep 17 00:00:00 2001 From: AngelFP Date: Tue, 16 Jun 2020 14:52:27 +0200 Subject: [PATCH 14/16] Make 'a' dimmensionless and avoid unit conversion --- .../data_handling/derived_field_definitions.py | 2 +- visualpic/data_handling/unit_converters.py | 15 +++++++-------- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/visualpic/data_handling/derived_field_definitions.py b/visualpic/data_handling/derived_field_definitions.py index fec3741..6a3cb9e 100644 --- a/visualpic/data_handling/derived_field_definitions.py +++ b/visualpic/data_handling/derived_field_definitions.py @@ -149,7 +149,7 @@ def calculate_norm_vector_pot(data_list, sim_geometry, sim_params): norm_vector_pot = {'name': 'a', - 'units': 'm_e*c^2/e', + 'units': '', 'requirements': {'1d': ['Ez'], '2dcartesian': ['Ez', 'Ex'], '3dcartesian': ['Ez', 'Ex', 'Ey'], diff --git a/visualpic/data_handling/unit_converters.py b/visualpic/data_handling/unit_converters.py index 3991b6b..4863794 100644 --- a/visualpic/data_handling/unit_converters.py +++ b/visualpic/data_handling/unit_converters.py @@ -7,6 +7,7 @@ License: GNU GPL-3.0. """ +import warnings import numpy as np import scipy.constants as ct @@ -70,11 +71,6 @@ 'MV': 1e-6} -norm_potential_conversion = {'V': ct.m_e * ct.c**2 / ct.e, - 'kV': ct.m_e * ct.c**2 / ct.e * 1e-3, - 'MV': ct.m_e * ct.c**2 / ct.e * 1e-6} - - class UnitConverter(): def __init__(self): """ Initialize unit converter. """ @@ -91,16 +87,19 @@ def __init__(self): 'T/m': bfieldgradient_conversion, 'm_e*c': momentum_conversion, 'W/m^2': intensity_conversion, - 'V': potential_conversion, - 'm_e*c^2/e': norm_potential_conversion} + 'V': potential_conversion} self.si_units = list(self.conversion_factors.keys()) def convert_field_units(self, field_data, field_md, target_field_units=None, target_axes_units=None, axes_to_convert=None, target_time_units=None): - + # Dimmensionless fields will not be converted. convert_field = target_field_units is not None + if convert_field and field_md['field']['units'] == '': + convert_field = False + warnings.warn('Field is dimmensionless. ' + 'No field unit conversion will be performed.') convert_axes = target_axes_units is not None convert_time = target_time_units is not None From 7e17779723c19911d3d1e5b91a86d4321dfcd1c5 Mon Sep 17 00:00:00 2001 From: AngelFP Date: Tue, 16 Jun 2020 15:02:27 +0200 Subject: [PATCH 15/16] Revert a (now unnecessary) change. --- visualpic/data_handling/unit_converters.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/visualpic/data_handling/unit_converters.py b/visualpic/data_handling/unit_converters.py index 4863794..a3c473d 100644 --- a/visualpic/data_handling/unit_converters.py +++ b/visualpic/data_handling/unit_converters.py @@ -123,9 +123,9 @@ def convert_field_units(self, field_data, field_md, axes_to_convert, convert_time) # convert field data to desired units - field_units = field_md['field']['units'] if convert_field and (target_field_units != 'SI' and - target_field_units != field_units): + target_field_units not in self.si_units): + field_units = field_md['field']['units'] field_data = self.convert_data(field_data, field_units, target_field_units) field_md['field']['units'] = target_field_units From a6b3d368cbbafeddcfe910b6a995958926ae6c28 Mon Sep 17 00:00:00 2001 From: AngelFP Date: Tue, 16 Jun 2020 15:05:06 +0200 Subject: [PATCH 16/16] Fix comments --- visualpic/data_handling/unit_converters.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/visualpic/data_handling/unit_converters.py b/visualpic/data_handling/unit_converters.py index a3c473d..3ead119 100644 --- a/visualpic/data_handling/unit_converters.py +++ b/visualpic/data_handling/unit_converters.py @@ -75,8 +75,7 @@ class UnitConverter(): def __init__(self): """ Initialize unit converter. """ # For convenience, due to their common use, the momentum units 'm_e*c' - # and the normalized potential units 'm_e*c^2/e' are also considered SI - # units. + # are also considered SI units. self.conversion_factors = {'m': length_conversion, 's': time_conversion, 'rad': angle_conversion, @@ -94,8 +93,8 @@ def __init__(self): def convert_field_units(self, field_data, field_md, target_field_units=None, target_axes_units=None, axes_to_convert=None, target_time_units=None): - # Dimmensionless fields will not be converted. convert_field = target_field_units is not None + # Dimmensionless fields will not be converted. if convert_field and field_md['field']['units'] == '': convert_field = False warnings.warn('Field is dimmensionless. '