Implementing HiPACE unit conversion definitions and some fixes

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,11 +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 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': '',
                    'requirements': {'1d': ['Ez'],
                                     '2dcartesian': ['Ez', 'Ex'],
                                     '3dcartesian': ['Ez', 'Ex', 'Ey'],

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
@@ -65,8 +66,16 @@
 intensity_conversion = {'W/cm^2': 1e-4}
 
 
+potential_conversion = {'m_e*c^2/e': 1 / ((ct.m_e * ct.c**2) / ct.e),
+                        'kV': 1e-3,
+                        'MV': 1e-6}
+
+
 class UnitConverter():
     def __init__(self):
+        """ Initialize unit converter. """
+        # For convenience, due to their common use, the momentum units 'm_e*c'
+        # are also considered SI units.
         self.conversion_factors = {'m': length_conversion,
                                    's': time_conversion,
                                    'rad': angle_conversion,
@@ -76,15 +85,20 @@ 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())
 
     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):
-
         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. '
+                          'No field unit conversion will be performed.')
         convert_axes = target_axes_units is not None
         convert_time = target_time_units is not None
 
@@ -243,8 +257,10 @@ 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']
+            }
         else:
             self.osiris_unit_conversion = None
         super().__init__()
@@ -272,7 +288,36 @@ 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.')

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":
@@ -180,7 +181,11 @@ 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])
+        attr_path = '/'
+        # In older Osiris versions the units are in field_path.
+        if "UNITS" in file[field_path].attrs:
+            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"""
@@ -197,26 +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 """
+        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.attrs['XMIN'][0],
-                                              file.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.attrs['XMIN'][1],
-                                                  file.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.attrs['XMIN'][2],
-                                                  file.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
 
@@ -264,6 +274,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()
@@ -274,9 +286,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 ''
 
@@ -288,17 +300,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)
@@ -308,7 +320,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
 
 

visualpic/data_reading/particle_readers.py

@@ -67,17 +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':
-                metadata['units'] = self._numpy_bytes_to_string(file_handle[
-                    self.name_relations[component]].attrs['UNITS'][0])
+                osiris_name = self.name_relations[component]
+                # In new Osiris versions, the units are in a list in '/'.
+                if 'QUANTS' in file_handle.attrs:
+                    units_path = '/'
+                    quantlist = [self._numpy_bytes_to_string(q)
+                                 for q in file_handle.attrs['QUANTS']]
+                    idx = quantlist.index(osiris_name)
+                else:
+                    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])
+            # 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 = '/'
+            sim_data = file_handle[simdata_path]
             metadata['grid'] = {}
-            metadata['grid']['resolution'] = file_handle.attrs['NX']
-            max_range = file_handle.attrs['XMAX']
-            min_range = file_handle.attrs['XMIN']
+            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):
@@ -112,26 +130,21 @@ def _read_component_data(self, file_path, species, component):
                 a = data[:, 0]
                 b = data[:, 1]
                 data = 1/2*(a+b)*(a+b+1)+b
-            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):
         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']
@@ -141,7 +154,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