From 12991c7f41451bb9440beadba6ed89ff7560485c Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 21 Apr 2020 07:33:25 -0700 Subject: [PATCH] Start implementing slice absolute position --- .../data_reader/field_reader.py | 35 +++++++++++++++---- openpmd_viewer/openpmd_timeseries/main.py | 29 ++++++++++----- .../openpmd_timeseries/utilities.py | 13 ++++++- 3 files changed, 62 insertions(+), 15 deletions(-) diff --git a/openpmd_viewer/openpmd_timeseries/data_reader/field_reader.py b/openpmd_viewer/openpmd_timeseries/data_reader/field_reader.py index ceb005fa..d09e7fe0 100644 --- a/openpmd_viewer/openpmd_timeseries/data_reader/field_reader.py +++ b/openpmd_viewer/openpmd_timeseries/data_reader/field_reader.py @@ -15,7 +15,7 @@ from ..utilities import construct_3d_from_circ def read_field_cartesian( filename, field_path, axis_labels, - slice_relative_position, slice_across ): + slice_relative_position, slice_absolute_position, slice_across ): """ Extract a given field from an HDF5 file in the openPMD format, when the geometry is cartesian (1d, 2d or 3d). @@ -47,6 +47,10 @@ def read_field_cartesian( filename, field_path, axis_labels, 0 : middle of the simulation box 1 : upper edge of the simulation box + slice_absolute_position : list of float, or None + Position(s), in meters, where to slice the data, along the + directions in `slice_across`s + Returns ------- A tuple with @@ -74,8 +78,16 @@ def read_field_cartesian( filename, field_path, axis_labels, list_slicing_index.append(slicing_index) # Number of cells along the slicing direction n_cells = shape[ slicing_index ] - # Index of the slice (prevent stepping out of the array) - i_cell = int( 0.5 * (slice_relative_position[count] + 1.) * n_cells ) + # Index of the slice + if slice_relative_position is not None: + i_cell = int( 0.5 * (slice_relative_position[count] + 1.) * n_cells ) + elif slice_absolute_position is not None: + i_cell = int( + (slice_absolute_position[count]-global_offset[slicing_index])/grid_spacing[slicing_index]) ) + else: + raise ValueError( 'slice_relative_position and ' + 'slice_absolute_position cannot both be None.') + # Prevent stepping out of the array i_cell = max( i_cell, 0 ) i_cell = min( i_cell, n_cells - 1) list_i_cell.append(i_cell) @@ -109,7 +121,7 @@ def read_field_cartesian( filename, field_path, axis_labels, def read_field_circ( filename, field_path, slice_relative_position, - slice_across, m=0, theta=0. ): + slice_absolute_position, slice_across, m=0, theta=0. ): """ Extract a given field from an HDF5 file in the openPMD format, when the geometry is thetaMode @@ -144,6 +156,10 @@ def read_field_circ( filename, field_path, slice_relative_position, 0 : middle of the simulation box 1 : upper edge of the simulation box + slice_absolute_position : list of float, or None + Position(s), in meters, where to slice the data, along the + directions in `slice_across`s + Returns ------- A tuple with @@ -231,8 +247,15 @@ def read_field_circ( filename, field_path, slice_relative_position, coord_array = getattr( info, slice_across_item ) # Number of cells along the slicing direction n_cells = len(coord_array) - # Index of the slice (prevent stepping out of the array) - i_cell = int( 0.5 * (slice_relative_position[count] + 1.) * n_cells ) + # Index of the slice + if slice_relative_position is not None: + i_cell = int( 0.5 * (slice_relative_position[count] + 1.) * n_cells ) + elif slice_absolute_position: + i_cell = np.argmin( abs(coord_array-slice_absolute_position[count]) ) + else: + raise ValueError( 'slice_relative_position and ' + 'slice_absolute_position cannot both be None.') + # Present stepping out of the array i_cell = max( i_cell, 0 ) i_cell = min( i_cell, n_cells - 1) F_total = np.take( F_total, [i_cell], axis=slicing_index ) diff --git a/openpmd_viewer/openpmd_timeseries/main.py b/openpmd_viewer/openpmd_timeseries/main.py index ca1bd747..21db0902 100644 --- a/openpmd_viewer/openpmd_timeseries/main.py +++ b/openpmd_viewer/openpmd_timeseries/main.py @@ -409,8 +409,16 @@ def get_field(self, field=None, coord=None, t=None, iteration=None, -1 : lower edge of the simulation box 0 : middle of the simulation box 1 : upper edge of the simulation box - Default: None, which results in slicing at 0 in all direction - of `slice_across`. + Default: None, which results in slicing at the middle of the box, + in each direction of `slice_across + (unless `slice_absolute_position` is set) + + slice_absolute_position : float or list of float, optional + Position(s), in meters, where to slice the data, along the + directions in `slice_across` + Default: None, which results in slicing at the middle of the box, + in each direction of `slice_across` + (unless `slice_relative_position` is set) plot : bool, optional Whether to plot the requested quantity @@ -442,8 +450,9 @@ def get_field(self, field=None, coord=None, t=None, iteration=None, "The available fields are: \n - %s\nPlease set the `field` " "argument accordingly." % field_list) # Check slicing - slice_across, slice_relative_position = \ - sanitize_slicing(slice_across, slice_relative_position) + slice_across, slice_relative_position, slice_absolute_position = \ + sanitize_slicing( slice_across, slice_relative_position, + slice_absolute_position ) if slice_across is not None: # Check that the elements are valid axis_labels = self.fields_metadata[field]['axis_labels'] @@ -495,21 +504,25 @@ def get_field(self, field=None, coord=None, t=None, iteration=None, # - For cartesian if geometry in ["1dcartesian", "2dcartesian", "3dcartesian"]: F, info = read_field_cartesian( filename, field_path, - axis_labels, slice_relative_position, slice_across) + axis_labels, slice_relative_position, slice_absolute_position, + slice_across) # - For thetaMode elif geometry == "thetaMode": if (coord in ['x', 'y']) and \ (self.fields_metadata[field]['type'] == 'vector'): # For Cartesian components, combine r and t components Fr, info = read_field_circ(filename, field + '/r', - slice_relative_position, slice_across, m, theta) + slice_relative_position, slice_absolute_position, + slice_across, m, theta) Ft, info = read_field_circ(filename, field + '/t', - slice_relative_position, slice_across, m, theta) + slice_relative_position, slice_absolute_position, + slice_across, m, theta) F = combine_cylindrical_components(Fr, Ft, theta, coord, info) else: # For cylindrical or scalar components, no special treatment F, info = read_field_circ(filename, field_path, - slice_relative_position, slice_across, m, theta) + slice_relative_position, slice_absolute_position, + slice_across, m, theta) # Plot the resulting field # Deactivate plotting when there is no slice selection diff --git a/openpmd_viewer/openpmd_timeseries/utilities.py b/openpmd_viewer/openpmd_timeseries/utilities.py index 6ce8cabd..6538b17e 100644 --- a/openpmd_viewer/openpmd_timeseries/utilities.py +++ b/openpmd_viewer/openpmd_timeseries/utilities.py @@ -16,7 +16,8 @@ from .data_reader.particle_reader import read_species_data from .numba_wrapper import jit -def sanitize_slicing(slice_across, slice_relative_position): +def sanitize_slicing(slice_across, slice_relative_position, + slice_absolute_position): """ Return standardized format for `slice_across` and `slice_relative_position`: - either `slice_across` and `slice_relative_position` are both `None` (no slicing) @@ -27,6 +28,8 @@ def sanitize_slicing(slice_across, slice_relative_position): ---------- slice_relative_position : float, or list of float, or None + slice_absolute_position : float, or list of float, or None + slice_across : str, or list of str, or None Direction(s) across which the data should be sliced """ @@ -34,6 +37,14 @@ def sanitize_slicing(slice_across, slice_relative_position): if slice_across is None or slice_across == []: return None, None + # Check that `slice_relative_position` and `slice_absolute_position` + # are not simultaneously provided: + if (slice_relative_position is not None) and \ + (slice_absolute_position is not None): + raise ValueError( + 'The arguments `slice_relative_position` and `slice_absolute_position` ' + 'cannot be used simultaneously.') + # Convert to lists if not isinstance(slice_across, list): slice_across = [slice_across]