From bea890b10fa75947633211cd1ad72a369d08c7a0 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Fri, 24 May 2024 21:57:03 +0200 Subject: [PATCH 01/18] particle count in run --- pyphare/pyphare/pharesee/hierarchy.py | 61 ++++++++++++++++++--------- pyphare/pyphare/pharesee/run.py | 4 ++ 2 files changed, 46 insertions(+), 19 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy.py index a7231bcac..3686cb2b7 100644 --- a/pyphare/pyphare/pharesee/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy.py @@ -227,19 +227,30 @@ class Patch: A patch represents a hyper-rectangular region of space """ - def __init__(self, patch_datas, patch_id=""): + def __init__(self, patch_datas, patch_id="", layout=None, attrs=None): """ :param patch_datas: a list of PatchData objects these are assumed to "belong" to the Patch so to share the same origin, mesh size and box. """ - pdata0 = list(patch_datas.values())[0] # 0 represents all others - self.layout = pdata0.layout - self.box = pdata0.layout.box - self.origin = pdata0.layout.origin - self.dl = pdata0.layout.dl - self.patch_datas = patch_datas - self.id = patch_id + if layout is not None: + self.layout = layout + self.box = layout.box + self.origin = layout.origin + self.dl = layout.dl + self.patch_datas = patch_datas + self.id = patch_id + + if len(patch_datas): + pdata0 = list(patch_datas.values())[0] # 0 represents all others + self.layout = pdata0.layout + self.box = pdata0.layout.box + self.origin = pdata0.layout.origin + self.dl = pdata0.layout.dl + self.patch_datas = patch_datas + self.id = patch_id + + self.attrs = attrs def __str__(self): return f"Patch: box( {self.box}), id({self.id})" @@ -1415,10 +1426,10 @@ def hierarchy_fromh5(h5_filename, time, hier, silent=True): if not silent: print("creating hierarchy from all times in file") times = list(data_file[h5_time_grp_key].keys()) - hier = hierarchy_fromh5(h5_filename, time=times[0], hier=hier) + hier = hierarchy_fromh5(h5_filename, time=times[0], hier=hier, silent=silent) if len(times) > 1: for t in times[1:]: - hierarchy_fromh5(h5_filename, t, hier) + hierarchy_fromh5(h5_filename, t, hier, silent=silent) return hier if create_from_one_time(time, hier): @@ -1441,16 +1452,21 @@ def hierarchy_fromh5(h5_filename, time, hier, silent=True): for pkey in h5_patch_lvl_grp.keys(): h5_patch_grp = h5_patch_lvl_grp[pkey] + patch_datas = {} + layout = make_layout(h5_patch_grp, lvl_cell_width, interp) if patch_has_datasets(h5_patch_grp): - patch_datas = {} - layout = make_layout(h5_patch_grp, lvl_cell_width, interp) add_to_patchdata(patch_datas, h5_patch_grp, basename, layout) - patches[ilvl].append( - Patch(patch_datas, h5_patch_grp.name.split("/")[-1]) + patches[ilvl].append( + Patch( + patch_datas, + h5_patch_grp.name.split("/")[-1], + layout=layout, + attrs={k: v for k, v in h5_patch_grp.attrs.items()}, ) + ) - patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl]) + patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl]) diag_hier = PatchHierarchy( patch_levels, domain_box, refinement_ratio, t, data_file @@ -1517,11 +1533,18 @@ def hierarchy_fromh5(h5_filename, time, hier, silent=True): for ipatch, pkey in enumerate(h5_time_grp[plvl_key].keys()): h5_patch_grp = h5_time_grp[plvl_key][pkey] + layout = make_layout(h5_patch_grp, lvl_cell_width, interp) + patch_datas = {} if patch_has_datasets(h5_patch_grp): - layout = make_layout(h5_patch_grp, lvl_cell_width, interp) - patch_datas = {} add_to_patchdata(patch_datas, h5_patch_grp, basename, layout) - lvl_patches.append(Patch(patch_datas)) + lvl_patches.append( + Patch( + patch_datas, + h5_patch_grp.name.split("/")[-1], + layout=layout, + attrs={k: v for k, v in h5_patch_grp.attrs.items()}, + ) + ) patch_levels[ilvl] = PatchLevel(ilvl, lvl_patches) @@ -1532,7 +1555,7 @@ def hierarchy_fromh5(h5_filename, time, hier, silent=True): if not silent: print("loading all times in existing hier") for time in data_file[h5_time_grp_key].keys(): - hier = hierarchy_fromh5(h5_filename, time, hier) + hier = hierarchy_fromh5(h5_filename, time, hier, silent=silent) return hier diff --git a/pyphare/pyphare/pharesee/run.py b/pyphare/pyphare/pharesee/run.py index 842027157..4fae4dedf 100644 --- a/pyphare/pyphare/pharesee/run.py +++ b/pyphare/pyphare/pharesee/run.py @@ -405,6 +405,10 @@ def filename(name): return hier return self._get_hierarchy(time, filename(pop_name), hier=hier) + def GetParticleCount(self, time): + c = self._get_hierarchy(time, "particle_count.h5") + return c + def GetMass(self, pop_name): list_of_qty = ["density", "flux", "domain", "levelGhost", "patchGhost"] list_of_mass = [] From cfc8511fa82659adae71de046b1ee2bcbe1acad5 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Sat, 25 May 2024 12:51:29 +0200 Subject: [PATCH 02/18] refactor the python hierarchy module --- pyphare/pyphare/pharesee/geometry.py | 3 +- pyphare/pyphare/pharesee/hierarchy.py | 1817 ----------------- .../pyphare/pharesee/hierarchy/__init__.py | 32 + pyphare/pyphare/pharesee/hierarchy/fromh5.py | 306 +++ pyphare/pyphare/pharesee/hierarchy/fromsim.py | 123 ++ .../pyphare/pharesee/hierarchy/hierarchy.py | 674 ++++++ .../pharesee/hierarchy/hierarchy_utils.py | 402 ++++ pyphare/pyphare/pharesee/hierarchy/patch.py | 66 + .../pyphare/pharesee/hierarchy/patchdata.py | 216 ++ .../pyphare/pharesee/hierarchy/patchlevel.py | 15 + pyphare/pyphare/pharesee/plotting.py | 2 +- pyphare/pyphare/pharesee/run.py | 9 +- .../pyphare_tests/test_pharesee/__init__.py | 7 +- .../test_pharesee/test_geometry.py | 9 +- .../test_pharesee/test_geometry_2d.py | 8 - .../field/coarsening/test_coarsen_field.py | 2 +- .../data/field/refine/test_refine_field.py | 2 +- tests/simulator/test_advance.py | 3 +- tests/simulator/test_diagnostics.py | 5 +- tests/simulator/test_initialization.py | 9 +- 20 files changed, 1860 insertions(+), 1850 deletions(-) delete mode 100644 pyphare/pyphare/pharesee/hierarchy.py create mode 100644 pyphare/pyphare/pharesee/hierarchy/__init__.py create mode 100644 pyphare/pyphare/pharesee/hierarchy/fromh5.py create mode 100644 pyphare/pyphare/pharesee/hierarchy/fromsim.py create mode 100644 pyphare/pyphare/pharesee/hierarchy/hierarchy.py create mode 100644 pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py create mode 100644 pyphare/pyphare/pharesee/hierarchy/patch.py create mode 100644 pyphare/pyphare/pharesee/hierarchy/patchdata.py create mode 100644 pyphare/pyphare/pharesee/hierarchy/patchlevel.py diff --git a/pyphare/pyphare/pharesee/geometry.py b/pyphare/pyphare/pharesee/geometry.py index 370f7da4f..a7cce9e1e 100644 --- a/pyphare/pyphare/pharesee/geometry.py +++ b/pyphare/pyphare/pharesee/geometry.py @@ -2,7 +2,8 @@ from ..core import box as boxm from pyphare.core.box import Box -from .hierarchy import FieldData, is_root_lvl +from .hierarchy.patchdata import FieldData +from .hierarchy.hierarchy_utils import is_root_lvl from pyphare.core.phare_utilities import listify, is_scalar diff --git a/pyphare/pyphare/pharesee/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy.py deleted file mode 100644 index 3686cb2b7..000000000 --- a/pyphare/pyphare/pharesee/hierarchy.py +++ /dev/null @@ -1,1817 +0,0 @@ -import os - -import numpy as np -import matplotlib.pyplot as plt - -from ..core import box as boxm -from ..core.box import Box -from ..core.gridlayout import GridLayout -from ..core.phare_utilities import deep_copy, refinement_ratio -from .particles import Particles -from ..core.phare_utilities import listify - - -class PatchData: - """ - base class for FieldData and ParticleData - this class just factors common geometrical properties - """ - - def __init__(self, layout, quantity): - """ - :param layout: a GridLayout representing the domain on which the data - is defined - :param quantity: ['field', 'particle'] - """ - self.quantity = quantity - self.box = layout.box - self.origin = layout.origin - self.layout = layout - - def __deepcopy__(self, memo): - no_copy_keys = ["dataset"] # do not copy these things - return deep_copy(self, memo, no_copy_keys) - - -class FieldData(PatchData): - """ - Concrete type of PatchData representing a physical quantity - defined on a grid. - """ - - @property - def x(self): - withGhosts = self.field_name != "tags" - if self._x is None: - self._x = self.layout.yeeCoordsFor( - self.field_name, - "x", - withGhosts=withGhosts, - centering=self.centerings[0], - ) - return self._x - - @property - def y(self): - withGhosts = self.field_name != "tags" - if self._y is None: - self._y = self.layout.yeeCoordsFor( - self.field_name, - "y", - withGhosts=withGhosts, - centering=self.centerings[1], - ) - return self._y - - @property - def z(self): - withGhosts = self.field_name != "tags" - if self._z is None: - self._z = self.layout.yeeCoordsFor( - self.field_name, - "z", - withGhosts=withGhosts, - centering=self.centerings[2], - ) - return self._z - - def primal_directions(self): - return self.size - self.ghost_box.shape - - def __str__(self): - return "FieldData: (box=({}, {}), key={})".format( - self.layout.box, self.layout.box.shape, self.field_name - ) - - def __repr__(self): - return self.__str__() - - def select(self, box): - """ - return view of internal data based on overlap of input box - returns a view +1 in size in primal directions - """ - assert isinstance(box, Box) and box.ndim == self.box.ndim - - gbox = self.ghost_box.copy() - gbox.upper += self.primal_directions() - - box = box.copy() - box.upper += self.primal_directions() - - overlap = box * gbox - if overlap is not None: - lower = self.layout.AMRToLocal(overlap.lower) - upper = self.layout.AMRToLocal(overlap.upper) - - if box.ndim == 1: - return self.dataset[lower[0] : upper[0] + 1] - if box.ndim == 2: - return self.dataset[lower[0] : upper[0] + 1, lower[1] : upper[1] + 1] - return np.array([]) - - def __getitem__(self, box): - return self.select(box) - - def __init__(self, layout, field_name, data, **kwargs): - """ - :param layout: A GridLayout representing the domain on which data is defined - :param field_name: the name of the field (e.g. "Bx") - :param data: the dataset from which data can be accessed - """ - super().__init__(layout, "field") - self._x = None - self._y = None - self._z = None - - self.layout = layout - self.field_name = field_name - self.name = field_name - self.dl = np.asarray(layout.dl) - self.ndim = layout.box.ndim - self.ghosts_nbr = np.zeros(self.ndim, dtype=int) - - if field_name in layout.centering["X"]: - directions = ["X", "Y", "Z"][: layout.box.ndim] # drop unused directions - self.centerings = [ - layout.qtyCentering(field_name, direction) for direction in directions - ] - elif "centering" in kwargs: - if isinstance(kwargs["centering"], list): - self.centerings = kwargs["centering"] - assert len(self.centerings) == self.ndim - else: - if self.ndim != 1: - raise ValueError( - "FieldData invalid dimenion for centering argument, expected list for dim > 1" - ) - self.centerings = [kwargs["centering"]] - else: - raise ValueError( - f"centering not specified and cannot be inferred from field name : {field_name}" - ) - - if self.field_name != "tags": - for i, centering in enumerate(self.centerings): - self.ghosts_nbr[i] = layout.nbrGhosts(layout.interp_order, centering) - - self.ghost_box = boxm.grow(layout.box, self.ghosts_nbr) - - self.size = np.copy(self.ghost_box.shape) - self.offset = np.zeros(self.ndim) - - for i, centering in enumerate(self.centerings): - if centering == "primal": - self.size[i] = self.ghost_box.shape[i] + 1 - else: - self.size[i] = self.ghost_box.shape[i] - self.offset[i] = 0.5 * self.dl[i] - - self.dataset = data - - def meshgrid(self, select=None): - def grid(): - if self.ndim == 1: - return [self.x] - if self.ndim == 2: - return np.meshgrid(self.x, self.y, indexing="ij") - return np.meshgrid(self.x, self.y, self.z, indexing="ij") - - mesh = grid() - if select is not None: - return tuple(g[select] for g in mesh) - return mesh - - -class ParticleData(PatchData): - """ - Concrete type of PatchData representing particles in a region - """ - - def __init__(self, layout, data, pop_name): - """ - :param layout: A GridLayout object representing the domain in which particles are - :param data: dataset containing particles - """ - super().__init__(layout, "particles") - self.dataset = data - self.pop_name = pop_name - self.name = pop_name - self.ndim = layout.box.ndim - - self.pop_name = pop_name - if layout.interp_order == 1: - self.ghosts_nbr = np.array([1] * layout.box.ndim) - elif layout.interp_order == 2 or layout.interp_order == 3: - self.ghosts_nbr = np.array([2] * layout.box.ndim) - else: - raise RuntimeError( - "invalid interpolation order {}".format(layout.interp_order) - ) - - self.ghost_box = boxm.grow(layout.box, self.ghosts_nbr) - assert (self.box.lower == self.ghost_box.lower + self.ghosts_nbr).all() - - def select(self, box): - return self.dataset[box] - - def __getitem__(self, box): - return self.select(box) - - def size(self): - return self.dataset.size() - - -class Patch: - """ - A patch represents a hyper-rectangular region of space - """ - - def __init__(self, patch_datas, patch_id="", layout=None, attrs=None): - """ - :param patch_datas: a list of PatchData objects - these are assumed to "belong" to the Patch so to - share the same origin, mesh size and box. - """ - if layout is not None: - self.layout = layout - self.box = layout.box - self.origin = layout.origin - self.dl = layout.dl - self.patch_datas = patch_datas - self.id = patch_id - - if len(patch_datas): - pdata0 = list(patch_datas.values())[0] # 0 represents all others - self.layout = pdata0.layout - self.box = pdata0.layout.box - self.origin = pdata0.layout.origin - self.dl = pdata0.layout.dl - self.patch_datas = patch_datas - self.id = patch_id - - self.attrs = attrs - - def __str__(self): - return f"Patch: box( {self.box}), id({self.id})" - - def __repr__(self): - return self.__str__() - - def copy(self): - """does not copy patchdatas.datasets (see class PatchData)""" - from copy import deepcopy - - return deepcopy(self) - - def __copy__(self): - return self.copy() - - def __call__(self, qty, **kwargs): - # take slice/slab of 1/2d array from 2/3d array - if "x" in kwargs and len(kwargs) == 1: - cut = kwargs["x"] - idim = 0 - elif "y" in kwargs and len(kwargs) == 1: - cut = kwargs["y"] - idim = 1 - else: - raise ValueError("need to specify either x or y cut coordinate") - pd = self.patch_datas[qty] - origin = pd.origin[idim] - idx = int((cut - origin) / pd.layout.dl[idim]) - nbrGhosts = pd.ghosts_nbr[idim] - if idim == 0: - return pd.dataset[idx + nbrGhosts, nbrGhosts:-nbrGhosts] - elif idim == 1: - return pd.dataset[nbrGhosts:-nbrGhosts, idx + nbrGhosts] - - -class PatchLevel: - """is a collection of patches""" - - def __init__(self, lvl_nbr, patches): - self.level_number = lvl_nbr - self.patches = patches - - def __iter__(self): - return self.patches.__iter__() - - def level_range(self): - name = list(self.patches[0].patch_datas.keys())[0] - return min([patch.patch_datas[name].x.min() for patch in self.patches]), max( - [patch.patch_datas[name].x.max() for patch in self.patches] - ) - - -def are_adjacent(lower, upper, atol=1e-6): - return np.abs(upper[0] - lower[-1]) < atol - - -def overlap_mask_1d(x, dl, level, qty): - """ - return the mask for x where x is overlaped by the qty patch datas - on the given level, assuming that this level is finer than the one of x - - :param x: 1d array containing the [x] position - :param dl: list containing the grid steps where x is defined - :param level: a given level associated to a hierarchy - :param qty: ['Bx', 'By', 'Bz', 'Ex', 'Ey', 'Ez', 'Fx', 'Fy', 'Fz', 'Vx', 'Vy', 'Vz', 'rho'] - """ - - is_overlaped = np.ones(x.shape[0], dtype=bool) * False - - for patch in level.patches: - pdata = patch.patch_datas[qty] - ghosts_nbr = pdata.ghosts_nbr - - fine_x = pdata.x[ghosts_nbr[0] - 1 : -ghosts_nbr[0] + 1] - - fine_dl = pdata.dl - local_dl = dl - - if fine_dl[0] < local_dl[0]: - xmin, xmax = fine_x.min(), fine_x.max() - - overlaped_idx = np.where((x > xmin) & (x < xmax))[0] - - is_overlaped[overlaped_idx] = True - - else: - raise ValueError("level needs to have finer grid resolution than that of x") - - return is_overlaped - - -def overlap_mask_2d(x, y, dl, level, qty): - """ - return the mask for x & y where ix & y are overlaped by the qty patch datas - on the given level, assuming that this level is finer than the one of x & y - important note : this mask is flatten - - :param x: 1d array containing the [x] position - :param y: 1d array containing the [y] position - :param dl: list containing the grid steps where x and y are defined - :param level: a given level associated to a hierarchy - :param qty: ['Bx', 'By', 'Bz', 'Ex', 'Ey', 'Ez', 'Fx', 'Fy', 'Fz', 'Vx', 'Vy', 'Vz', 'rho'] - """ - - is_overlaped = np.ones([x.shape[0] * y.shape[0]], dtype=bool) * False - - for patch in level.patches: - pdata = patch.patch_datas[qty] - ghosts_nbr = pdata.ghosts_nbr - - fine_x = pdata.x[ghosts_nbr[0] - 1 : -ghosts_nbr[0] + 1] - fine_y = pdata.y[ghosts_nbr[1] - 1 : -ghosts_nbr[1] + 1] - - fine_dl = pdata.dl - local_dl = dl - - if (fine_dl[0] < local_dl[0]) and (fine_dl[1] < local_dl[1]): - xmin, xmax = fine_x.min(), fine_x.max() - ymin, ymax = fine_y.min(), fine_y.max() - - xv, yv = np.meshgrid(x, y, indexing="ij") - xf = xv.flatten() - yf = yv.flatten() - - overlaped_idx = np.where( - (xf > xmin) & (xf < xmax) & (yf > ymin) & (yf < ymax) - )[0] - - is_overlaped[overlaped_idx] = True - - else: - raise ValueError( - "level needs to have finer grid resolution than that of x or y" - ) - - return is_overlaped - - -def flat_finest_field(hierarchy, qty, time=None): - """ - returns 2 flattened arrays containing the data (with shape [Npoints]) - and the coordinates (with shape [Npoints, Ndim]) for the given - hierarchy of qty. - - :param hierarchy: the hierarchy where qty is defined - :param qty: the field (eg "Bx") that we want - """ - - dim = hierarchy.ndim - - if dim == 1: - return flat_finest_field_1d(hierarchy, qty, time) - elif dim == 2: - return flat_finest_field_2d(hierarchy, qty, time) - elif dim == 3: - raise RuntimeError("Not implemented") - # return flat_finest_field_3d(hierarchy, qty, time) - else: - raise ValueError("the dim of a hierarchy should be 1, 2 or 3") - - -def flat_finest_field_1d(hierarchy, qty, time=None): - lvl = hierarchy.levels(time) - - for ilvl in range(hierarchy.finest_level(time) + 1)[::-1]: - patches = lvl[ilvl].patches - - for ip, patch in enumerate(patches): - pdata = patch.patch_datas[qty] - - # all but 1 ghost nodes are removed in order to limit - # the overlapping, but to keep enough point to avoid - # any extrapolation for the interpolator - needed_points = pdata.ghosts_nbr - 1 - - # data = pdata.dataset[patch.box] # TODO : once PR 551 will be merged... - data = pdata.dataset[needed_points[0] : -needed_points[0]] - x = pdata.x[needed_points[0] : -needed_points[0]] - - if ilvl == hierarchy.finest_level(time): - if ip == 0: - final_data = data - final_x = x - else: - final_data = np.concatenate((final_data, data)) - final_x = np.concatenate((final_x, x)) - - else: - is_overlaped = overlap_mask_1d( - x, pdata.dl, hierarchy.level(ilvl + 1, time), qty - ) - - finest_data = data[~is_overlaped] - finest_x = x[~is_overlaped] - - final_data = np.concatenate((final_data, finest_data)) - final_x = np.concatenate((final_x, finest_x)) - - return final_data, final_x - - -def flat_finest_field_2d(hierarchy, qty, time=None): - lvl = hierarchy.levels(time) - - for ilvl in range(hierarchy.finest_level(time) + 1)[::-1]: - patches = lvl[ilvl].patches - - for ip, patch in enumerate(patches): - pdata = patch.patch_datas[qty] - - # all but 1 ghost nodes are removed in order to limit - # the overlapping, but to keep enough point to avoid - # any extrapolation for the interpolator - needed_points = pdata.ghosts_nbr - 1 - - # data = pdata.dataset[patch.box] # TODO : once PR 551 will be merged... - data = pdata.dataset[ - needed_points[0] : -needed_points[0], - needed_points[1] : -needed_points[1], - ] - x = pdata.x[needed_points[0] : -needed_points[0]] - y = pdata.y[needed_points[1] : -needed_points[1]] - - xv, yv = np.meshgrid(x, y, indexing="ij") - - data_f = data.flatten() - xv_f = xv.flatten() - yv_f = yv.flatten() - - if ilvl == hierarchy.finest_level(time): - if ip == 0: - final_data = data_f - tmp_x = xv_f - tmp_y = yv_f - else: - final_data = np.concatenate((final_data, data_f)) - tmp_x = np.concatenate((tmp_x, xv_f)) - tmp_y = np.concatenate((tmp_y, yv_f)) - - else: - is_overlaped = overlap_mask_2d( - x, y, pdata.dl, hierarchy.level(ilvl + 1, time), qty - ) - - finest_data = data_f[~is_overlaped] - finest_x = xv_f[~is_overlaped] - finest_y = yv_f[~is_overlaped] - - final_data = np.concatenate((final_data, finest_data)) - tmp_x = np.concatenate((tmp_x, finest_x)) - tmp_y = np.concatenate((tmp_y, finest_y)) - - final_xy = np.stack((tmp_x, tmp_y), axis=1) - - return final_data, final_xy - - -def finest_part_data(hierarchy, time=None): - """ - returns a dict {popname : Particles} - Particles contained in the dict are those from - the finest patches available at a given location - """ - from copy import deepcopy - - from .particles import remove - - # we are going to return a dict {popname : Particles} - # we prepare it with population names - aPatch = hierarchy.level(0, time=time).patches[0] - particles = {popname: None for popname in aPatch.patch_datas.keys()} - - # our strategy is to explore the hierarchy from the finest - # level to the coarsest. at Each level we keep only particles - # that are in cells that are not overlaped by finer boxes - - # this dict keeps boxes for patches at each level - # each level will thus need this dict to see next finer boxes - lvlPatchBoxes = {ilvl: [] for ilvl in range(hierarchy.finest_level(time) + 1)} - - for ilvl in range(hierarchy.finest_level(time) + 1)[::-1]: - plvl = hierarchy.level(ilvl, time=time) - for ip, patch in enumerate(plvl.patches): - lvlPatchBoxes[ilvl].append(patch.box) - for popname, pdata in patch.patch_datas.items(): - # if we're at the finest level - # we need to keep all particles - if ilvl == hierarchy.finest_level(time): - if particles[popname] is None: - particles[popname] = deepcopy(pdata.dataset) - else: - particles[popname].add(deepcopy(pdata.dataset)) - - # if there is a finer level - # we need to keep only those of the current patch - # that are not in cells overlaped by finer patch boxes - else: - icells = pdata.dataset.iCells - parts = deepcopy(pdata.dataset) - create = True - for finerBox in lvlPatchBoxes[ilvl + 1]: - coarseFinerBox = boxm.coarsen(finerBox, refinement_ratio) - within = np.where( - (icells >= coarseFinerBox.lower[0]) - & (icells <= coarseFinerBox.upper[0]) - )[0] - if create: - toRemove = within - create = False - else: - toRemove = np.concatenate((toRemove, within)) - - if toRemove.size != 0: - parts = remove(parts, toRemove) - if parts is not None: - particles[popname].add(parts) - return particles - - -class PatchHierarchy(object): - """is a collection of patch levels""" - - def __init__( - self, - patch_levels, - domain_box, - refinement_ratio=2, - time=0.0, - data_files=None, - **kwargs, - ): - self.patch_levels = patch_levels - self.ndim = len(domain_box.lower) - self.time_hier = {} - self.time_hier.update({self.format_timestamp(time): patch_levels}) - - self.domain_box = domain_box - self.refinement_ratio = refinement_ratio - - self.data_files = {} - self._sim = None - - if data_files is not None: - self.data_files.update(data_files) - - if len(self.quantities()) > 1: - for qty in self.quantities(): - if qty in self.__dict__: - continue - first = True - for time, levels in self.time_hier.items(): - new_lvls = {} - for ilvl, level in levels.items(): - patches = [] - for patch in level.patches: - patches += [Patch({qty: patch.patch_datas[qty]}, patch.id)] - new_lvls[ilvl] = PatchLevel(ilvl, patches) - if first: - self.__dict__[qty] = PatchHierarchy( - new_lvls, self.domain_box, time=time - ) - first = False - else: - self.qty.time_hier[time] = new_lvls # pylint: disable=E1101 - - def __getitem__(self, qty): - return self.__dict__[qty] - - @property - def sim(self): - if self._sim: - return self._sim - - if "py_attrs" not in self.data_files: - raise ValueError("Simulation is not available for deserialization") - - from ..pharein.simulation import deserialize - - try: - self._sim = deserialize( - self.data_files["py_attrs"].attrs["serialized_simulation"] - ) - except Exception as e: - raise RuntimeError(f"Failed to deserialize simulation from data file : {e}") - return self._sim - - def __call__(self, qty=None, **kwargs): - # take slice/slab of 1/2d array from 2/3d array - def cuts(c, coord): - return c > coord.min() and c < coord.max() - - class Extractor: - def __init__(self): - self.exclusions = [] - - def extract(self, coord, data): - mask = coord == coord - for exclusion in self.exclusions: - idx = np.where( - (coord > exclusion[0] - 1e-6) & (coord < exclusion[1] + 1e-6) - )[0] - mask[idx] = False - - self.exclusions += [(coord.min(), coord.max())] - return coord[mask], data[mask] - - def domain_coords(patch, qty): - pd = patch.patch_datas[qty] - nbrGhosts = pd.ghosts_nbr[0] - return pd.x[nbrGhosts:-nbrGhosts], pd.y[nbrGhosts:-nbrGhosts] - - if len(kwargs) < 1 or len(kwargs) > 3: - raise ValueError("Error - must provide coordinates") - if qty == None: - if len(self.quantities()) == 1: - qty = self.quantities()[0] - else: - raise ValueError( - "The PatchHierarchy has several quantities but none is specified" - ) - - if len(kwargs) == 1: # 1D cut - if "x" in kwargs: - c = kwargs["x"] - slice_dim = 1 - cst_dim = 0 - else: - c = kwargs["y"] - slice_dim = 0 - cst_dim = 1 - - extractor = Extractor() - datas = [] - coords = [] - ilvls = list(self.levels().keys())[::-1] - - for ilvl in ilvls: - lvl = self.patch_levels[ilvl] - for patch in lvl.patches: - slice_coord = domain_coords(patch, qty)[slice_dim] - cst_coord = domain_coords(patch, qty)[cst_dim] - - if cuts(c, cst_coord): - data = patch(qty, **kwargs) - coord_keep, data_keep = extractor.extract(slice_coord, data) - datas += [data_keep] - coords += [coord_keep] - - cut = np.concatenate(datas) - coords = np.concatenate(coords) - ic = np.argsort(coords) - coords = coords[ic] - cut = cut[ic] - return coords, cut - - def _default_time(self): - return self.times()[0] - - def finest_level(self, time=None): - if time is None: - time = self._default_time() - return max(list(self.levels(time=time).keys())) - - def levels(self, time=None): - if time is None: - time = self._default_time() - return self.time_hier[self.format_timestamp(time)] - - def level(self, level_number, time=None): - return self.levels(time)[level_number] - - def levelNbr(self, time=None): - if time is None: - time = self._default_time() - return len(self.levels(time).items()) - - def levelNbrs(self, time=None): - if time is None: - time = self._default_time() - return list(self.levels(time).keys()) - - def is_homogeneous(self): - """ - return True if all patches of all levels at all times - have the same patch data quantities - """ - qties = self._quantities() - it_is = True - for time, levels in self.time_hier.items(): - for ilvl, lvl in levels.items(): - for patch in lvl.patches: - it_is &= qties == list(patch.patch_datas.keys()) - return it_is - - def _quantities(self): - for ilvl, lvl in self.levels().items(): - if len(lvl.patches) > 0: - return list(self.level(ilvl).patches[0].patch_datas.keys()) - return [] - - def quantities(self): - if not self.is_homogeneous(): - raise RuntimeError("Error - hierarchy is not homogeneous") - return self._quantities() - - def global_min(self, qty, **kwargs): - time = kwargs.get("time", self._default_time()) - first = True - for ilvl, lvl in self.levels(time).items(): - for patch in lvl.patches: - pd = patch.patch_datas[qty] - if first: - m = pd.dataset[:].min() - first = False - else: - m = min(m, pd.dataset[:].min()) - - return m - - def global_max(self, qty, **kwargs): - time = kwargs.get("time", self._default_time()) - first = True - for ilvl, lvl in self.levels(time).items(): - for patch in lvl.patches: - pd = patch.patch_datas[qty] - if first: - m = pd.dataset[:].max() - first = False - else: - m = max(m, pd.dataset[:].max()) - - return m - - def refined_domain_box(self, level_number): - """ - returns the domain box refined for a given level number - """ - assert level_number >= 0 - return boxm.refine(self.domain_box, self.refinement_ratio**level_number) - - def format_timestamp(self, timestamp): - if isinstance(timestamp, str): - return timestamp - return "{:.10f}".format(timestamp) - - def level_domain_box(self, level_number): - if level_number == 0: - return self.domain_box - return self.refined_domain_box(level_number) - - def __str__(self): - s = "Hierarchy: \n" - for t, patch_levels in self.time_hier.items(): - for ilvl, lvl in patch_levels.items(): - s = s + "Level {}\n".format(ilvl) - for ip, patch in enumerate(lvl.patches): - for qty_name, pd in patch.patch_datas.items(): - pdstr = " P{ip} {type} {pdname} box is {box} and ghost box is {gbox}" - s = s + pdstr.format( - ip=ip, - type=type(pd.dataset), - pdname=qty_name, - box=patch.box, - gbox=pd.ghost_box, - ) - s = s + "\n" - return s - - def times(self): - return np.sort(np.asarray(list(self.time_hier.keys()))) - - def plot_patches(self, save=False): - fig, ax = plt.subplots(figsize=(10, 3)) - for ilvl, lvl in self.levels(0.0).items(): - lvl_offset = ilvl * 0.1 - for patch in lvl.patches: - dx = patch.dl[0] - x0 = patch.box.lower * dx - x1 = patch.box.upper * dx - xcells = np.arange(x0, x1 + dx, dx) - y = lvl_offset + np.zeros_like(xcells) - ax.plot(xcells, y, marker=".") - - if save: - fig.savefig("hierarchy.png") - - def box_to_Rectangle(self, box): - from matplotlib.patches import Rectangle - - return Rectangle(box.lower, *box.shape) - - def plot_2d_patches(self, ilvl, collections, **kwargs): - if isinstance(collections, list) and all( - [isinstance(el, Box) for el in collections] - ): - collections = [{"boxes": collections}] - - from matplotlib.collections import PatchCollection - - level_domain_box = self.level_domain_box(ilvl) - mi, ma = level_domain_box.lower.min(), level_domain_box.upper.max() - - fig, ax = kwargs.get("subplot", plt.subplots(figsize=(6, 6))) - - for collection in collections: - facecolor = collection.get("facecolor", "none") - edgecolor = collection.get("edgecolor", "purple") - alpha = collection.get("alpha", 1) - rects = [self.box_to_Rectangle(box) for box in collection["boxes"]] - - ax.add_collection( - PatchCollection( - rects, facecolor=facecolor, alpha=alpha, edgecolor=edgecolor - ) - ) - - if "title" in kwargs: - from textwrap import wrap - - xfigsize = int(fig.get_size_inches()[0] * 10) # 10 characters per inch - ax.set_title("\n".join(wrap(kwargs["title"], xfigsize))) - - major_ticks = np.arange(mi - 5, ma + 5 + 5, 5) - ax.set_xticks(major_ticks) - ax.set_yticks(major_ticks) - - minor_ticks = np.arange(mi - 5, ma + 5 + 5, 1) - ax.set_xticks(minor_ticks, minor=True) - ax.set_yticks(minor_ticks, minor=True) - - ax.grid(which="both") - - return fig - - def plot1d(self, **kwargs): - """ - plot - """ - usr_lvls = kwargs.get("levels", (0,)) - qty = kwargs.get("qty", None) - time = kwargs.get("time", self.times()[0]) - - if "ax" not in kwargs: - fig, ax = plt.subplots() - else: - ax = kwargs["ax"] - fig = ax.figure - for lvl_nbr, level in self.levels(time).items(): - if lvl_nbr not in usr_lvls: - continue - for ip, patch in enumerate(level.patches): - pdata_nbr = len(patch.patch_datas) - pdata_names = list(patch.patch_datas.keys()) - if qty is None and pdata_nbr != 1: - multiple = "multiple quantities in patch, " - err = ( - multiple - + "please specify a quantity in " - + " ".join(pdata_names) - ) - raise ValueError(err) - if qty is None: - qty = pdata_names[0] - - layout = patch.patch_datas[qty].layout - nbrGhosts = layout.nbrGhostFor(qty) - val = patch.patch_datas[qty][patch.box] - x = patch.patch_datas[qty].x[nbrGhosts[0] : -nbrGhosts[0]] - label = "L{level}P{patch}".format(level=lvl_nbr, patch=ip) - marker = kwargs.get("marker", "") - ls = kwargs.get("ls", "--") - color = kwargs.get("color", "k") - ax.plot(x, val, label=label, marker=marker, ls=ls, color=color) - - ax.set_title(kwargs.get("title", "")) - ax.set_xlabel(kwargs.get("xlabel", "x")) - ax.set_ylabel(kwargs.get("ylabel", qty)) - if "xlim" in kwargs: - ax.set_xlim(kwargs["xlim"]) - if "ylim" in kwargs: - ax.set_ylim(kwargs["ylim"]) - - if kwargs.get("legend", None) is not None: - ax.legend() - - if "filename" in kwargs: - fig.savefig(kwargs["filename"]) - - def plot2d(self, **kwargs): - from matplotlib.patches import Rectangle - from mpl_toolkits.axes_grid1 import make_axes_locatable - - time = kwargs.get("time", self._default_time()) - usr_lvls = kwargs.get("levels", self.levelNbrs(time)) - default_qty = None - if len(self.quantities()) == 1: - default_qty = self.quantities()[0] - qty = kwargs.get("qty", default_qty) - - if "ax" not in kwargs: - fig, ax = plt.subplots() - else: - ax = kwargs["ax"] - fig = ax.figure - - glob_min = self.global_min(qty) - glob_max = self.global_max(qty) - # assumes max 5 levels... - patchcolors = {ilvl: "k" for ilvl in usr_lvls} - patchcolors = kwargs.get("patchcolors", patchcolors) - if not isinstance(patchcolors, dict): - patchcolors = dict(zip(usr_lvls, patchcolors)) - - linewidths = {ilvl: 1 for ilvl in usr_lvls} - linewidths = kwargs.get("lw", linewidths) - if not isinstance(linewidths, dict): - linewidths = dict(zip(usr_lvls, linewidths)) - linestyles = {ilvl: "-" for ilvl in usr_lvls} - linestyles = kwargs.get("ls", linestyles) - if not isinstance(linestyles, dict): - linestyles = dict(zip(usr_lvls, linestyles)) - - for lvl_nbr, lvl in self.levels(time).items(): - if lvl_nbr not in usr_lvls: - continue - for patch in self.level(lvl_nbr, time).patches: - pdat = patch.patch_datas[qty] - data = pdat.dataset[:] - nbrGhosts = pdat.ghosts_nbr - x = pdat.x - y = pdat.y - - # if nbrGhosts is 0, we cannot do array[0,-0] - if np.all(nbrGhosts == np.zeros_like(nbrGhosts)): - x = np.copy(x) - y = np.copy(y) - else: - data = pdat[patch.box] - x = np.copy(x[nbrGhosts[0] : -nbrGhosts[0]]) - y = np.copy(y[nbrGhosts[1] : -nbrGhosts[1]]) - dx, dy = pdat.layout.dl - x -= dx * 0.5 - y -= dy * 0.5 - x = np.append(x, x[-1] + dx) - y = np.append(y, y[-1] + dy) - im = ax.pcolormesh( - x, - y, - data.T, - cmap=kwargs.get("cmap", "Spectral_r"), - vmin=kwargs.get("vmin", glob_min - 1e-6), - vmax=kwargs.get("vmax", glob_max + 1e-6), - ) - - if kwargs.get("plot_patches", False) is True: - r = Rectangle( - (patch.box.lower[0] * dx, patch.box.lower[1] * dy), - patch.box.shape[0] * dx, - patch.box.shape[1] * dy, - fc="none", - ec=patchcolors[lvl_nbr], - alpha=0.4, - lw=linewidths[lvl_nbr], - ls=linestyles[lvl_nbr], - ) - ax.add_patch(r) - - ax.set_aspect(kwargs.get("aspect", "equal")) - ax.set_title(kwargs.get("title", "")) - ax.set_xlabel(kwargs.get("xlabel", "x")) - ax.set_ylabel(kwargs.get("ylabel", "y")) - if "xlim" in kwargs: - ax.set_xlim(kwargs["xlim"]) - if "ylim" in kwargs: - ax.set_ylim(kwargs["ylim"]) - - divider = make_axes_locatable(ax) - cax = divider.append_axes("right", size="5%", pad=0.08) - plt.colorbar(im, ax=ax, cax=cax) - - if kwargs.get("legend", None) is not None: - ax.legend() - - if "filename" in kwargs: - fig.savefig(kwargs["filename"], dpi=kwargs.get("dpi", 200)) - - return fig, ax - - def plot(self, **kwargs): - if self.ndim == 1: - return self.plot1d(**kwargs) - elif self.ndim == 2: - return self.plot2d(**kwargs) - - def dist_plot(self, **kwargs): - """ - plot phase space of a particle hierarchy - """ - import copy - - from .plotting import dist_plot as dp - - usr_lvls = kwargs.get("levels", (0,)) - finest = kwargs.get("finest", False) - pops = kwargs.get("pop", []) - time = kwargs.get("time", self.times()[0]) - axis = kwargs.get("axis", ("Vx", "Vy")) - all_pops = list(self.level(0, time).patches[0].patch_datas.keys()) - - vmin = kwargs.get("vmin", -2) - vmax = kwargs.get("vmax", 2) - dv = kwargs.get("dv", 0.05) - vbins = vmin + dv * np.arange(int((vmax - vmin) / dv)) - - if finest: - final = finest_part_data(self) - if axis[0] == "x": - xbins = amr_grid(self, time) - bins = (xbins, vbins) - else: - bins = (vbins, vbins) - kwargs["bins"] = bins - - else: - final = {pop: None for pop in all_pops} - for lvl_nbr, level in self.levels(time).items(): - if lvl_nbr not in usr_lvls: - continue - for ip, patch in enumerate(level.patches): - if len(pops) == 0: - pops = list(patch.patch_datas.keys()) - - for pop in pops: - tmp = copy.copy(patch.patch_datas[pop].dataset) - - if final[pop] is None: - final[pop] = tmp - else: - final[pop].add(tmp) - - # select particles - if "select" in kwargs: - for pop, particles in final.items(): - final[pop] = kwargs["select"](particles) - - return final, dp(final, **kwargs) - - -def amr_grid(hierarchy, time): - """returns a non-uniform contiguous primal grid - associated to the given hierarchy - """ - lvlPatchBoxes = {ilvl: [] for ilvl in range(hierarchy.finest_level() + 1)} - finalCells = {ilvl: None for ilvl in range(hierarchy.finest_level() + 1)} - lvl = hierarchy.levels(time) - - for ilvl in range(hierarchy.finest_level(time) + 1)[::-1]: - sorted_patches = sorted(lvl[ilvl].patches, key=lambda p: p.layout.box.lower[0]) - - for ip, patch in enumerate(sorted_patches): - box = patch.layout.box - lvlPatchBoxes[ilvl].append(box) - - # we create a list of all cells in the current patch - # remember that if the box upper cell is, say = 40, - # it means that the upper node is the lower node of cell 41 - # so to get all primal nodes of a patch we need to include - # one past the upper cell. - # this said we do not want to include that last primal nodes - # all the time because that would be a duplicate with the lower - # node of the next patch. We only want to add it for the LAST - # (because sorted) patch. We also do not want to do it on levels - # other than L0 because the last primal node of the last patch - # of L_i is the first primal node of a L_{i-1} node, so including it - # would also mean adding a duplicate. - last = 1 if ilvl == 0 and ip == len(sorted_patches) - 1 else 0 - cells = np.arange(box.lower[0], box.upper[0] + 1 + last) - - # finest level has no next finer so we take all cells - if ilvl == hierarchy.finest_level(time): - if finalCells[ilvl] is None: - finalCells[ilvl] = cells - else: - finalCells[ilvl] = np.concatenate((finalCells[ilvl], cells)) - - else: - # on other levels - # we take only grids not overlaped by next finer - coarsenedNextFinerBoxes = [ - boxm.coarsen(b, refinement_ratio) for b in lvlPatchBoxes[ilvl + 1] - ] - for coarseBox in coarsenedNextFinerBoxes: - ccells = np.arange(coarseBox.lower[0], coarseBox.upper[0] + 1) - inter, icells, iccells = np.intersect1d( - cells, ccells, return_indices=True - ) - cells = np.delete(cells, icells) - if len(cells): - if finalCells[ilvl] is None: - finalCells[ilvl] = cells - else: - finalCells[ilvl] = np.unique( - np.concatenate((finalCells[ilvl], cells)) - ) - - # now we have all cells for each level we - # just need to compute the primal coordinates - # and concatenate in a single array - for ilvl in range(hierarchy.finest_level() + 1): - if ilvl == 0: - x = finalCells[ilvl] * hierarchy.level(ilvl).patches[0].layout.dl[0] - else: - xx = finalCells[ilvl] * hierarchy.level(ilvl).patches[0].layout.dl[0] - x = np.concatenate((x, xx)) - - return np.sort(x) - - -def is_root_lvl(patch_level): - return patch_level.level_number == 0 - - -field_qties = { - "EM_B_x": "Bx", - "EM_B_y": "By", - "EM_B_z": "Bz", - "EM_E_x": "Ex", - "EM_E_y": "Ey", - "EM_E_z": "Ez", - "flux_x": "Fx", - "flux_y": "Fy", - "flux_z": "Fz", - "bulkVelocity_x": "Vx", - "bulkVelocity_y": "Vy", - "bulkVelocity_z": "Vz", - "momentum_tensor_xx": "Mxx", - "momentum_tensor_xy": "Mxy", - "momentum_tensor_xz": "Mxz", - "momentum_tensor_yy": "Myy", - "momentum_tensor_yz": "Myz", - "momentum_tensor_zz": "Mzz", - "density": "rho", - "mass_density": "rho", - "tags": "tags", -} - - -particle_files_patterns = ("domain", "patchGhost", "levelGhost") - - -def is_particle_file(filename): - return any([pattern in filename for pattern in particle_files_patterns]) - - -def particle_dataset_name(basename): - """ - return "alpha_domain" from ions_pop_alpha_domain.h5 - """ - popname = basename.strip(".h5").split("_")[-2] - part_type = basename.strip(".h5").split("_")[-1] - dataset_name = popname + "_" + part_type - - return dataset_name - - -def is_pop_fluid_file(basename): - return (is_particle_file(basename) is False) and "pop" in basename - - -def make_layout(h5_patch_grp, cell_width, interp_order): - origin = h5_patch_grp.attrs["origin"] - upper = h5_patch_grp.attrs["upper"] - lower = h5_patch_grp.attrs["lower"] - return GridLayout(Box(lower, upper), origin, cell_width, interp_order=interp_order) - - -def create_from_all_times(time, hier): - return time is None and hier is None - - -def create_from_one_time(time, hier): - return time is not None and hier is None - - -def load_all_times(time, hier): - return time is None and hier is not None - - -def load_one_time(time, hier): - return time is not None and hier is not None - - -def compute_hier_from_(h, compute): - """ - returns a hierarchy resulting from calling 'compute' - on each patch of the given hierarchy 'h' - - compute is a function taking a Patch and returning - a list of dicts with the following keys: - - "data": ndarray containing the data - "name": str, name of the data living on that patch, must be unique - "centering": str, ["dual", "primal"] - - caveat: routine only works in 1D so far. - """ - assert len(h.time_hier) == 1 # only single time hierarchies now - patch_levels = {} - for ilvl, lvl in h.patch_levels.items(): - patches = {} - for ip, patch in enumerate(lvl.patches): - new_patch_datas = {} - layout = patch.layout - datas = compute(patch) - for data in datas: - pd = FieldData( - layout, data["name"], data["data"], centering=data["centering"] - ) - new_patch_datas[data["name"]] = pd - if ilvl not in patches: - patches[ilvl] = [] - patches[ilvl].append(Patch(new_patch_datas, patch.id)) - - patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl]) - - t = list(h.time_hier.keys())[0] - return PatchHierarchy(patch_levels, h.domain_box, refinement_ratio, time=t) - - -def are_compatible_hierarchies(hierarchies): - return True - - -def extract_patchdatas(hierarchies, ilvl, ipatch): - """ - returns a dict {patchdata_name:patchdata} from a list of hierarchies for patch ipatch at level ilvl - """ - patches = [h.patch_levels[ilvl].patches[ipatch] for h in hierarchies] - patch_datas = { - pdname: pd for p in patches for pdname, pd in list(p.patch_datas.items()) - } - return patch_datas - - -def new_patchdatas_from(compute, patchdatas, layout, id, **kwargs): - new_patch_datas = {} - datas = compute(patchdatas, id=id, **kwargs) - for data in datas: - pd = FieldData(layout, data["name"], data["data"], centering=data["centering"]) - new_patch_datas[data["name"]] = pd - return new_patch_datas - - -def new_patches_from(compute, hierarchies, ilvl, **kwargs): - reference_hier = hierarchies[0] - new_patches = [] - patch_nbr = len(reference_hier.patch_levels[ilvl].patches) - for ip in range(patch_nbr): - current_patch = reference_hier.patch_levels[ilvl].patches[ip] - layout = current_patch.layout - patch_datas = extract_patchdatas(hierarchies, ilvl, ip) - new_patch_datas = new_patchdatas_from( - compute, patch_datas, layout, id=current_patch.id, **kwargs - ) - new_patches.append(Patch(new_patch_datas, current_patch.id)) - return new_patches - - -def compute_hier_from(compute, hierarchies, **kwargs): - """return a new hierarchy using the callback 'compute' on all patchdatas - of the given hierarchies - """ - if not are_compatible_hierarchies(hierarchies): - raise RuntimeError("hierarchies are not compatible") - - hierarchies = listify(hierarchies) - reference_hier = hierarchies[0] - domain_box = reference_hier.domain_box - patch_levels = {} - for ilvl in range(reference_hier.levelNbr()): - patch_levels[ilvl] = PatchLevel( - ilvl, new_patches_from(compute, hierarchies, ilvl, **kwargs) - ) - - assert len(reference_hier.time_hier) == 1 # only single time hierarchies now - t = list(reference_hier.time_hier.keys())[0] - return PatchHierarchy(patch_levels, domain_box, refinement_ratio, time=t) - - -def pop_name(basename): - return basename.strip(".h5").split("_")[2] - - -def add_to_patchdata(patch_datas, h5_patch_grp, basename, layout): - """ - adds data in the h5_patch_grp in the given PatchData dict - returns True if valid h5 patch found - """ - - if is_particle_file(basename): - v = np.asarray(h5_patch_grp["v"]) - s = v.size - v = v[:].reshape(int(s / 3), 3) - nbrParts = v.shape[0] - dl = np.zeros((nbrParts, layout.ndim)) - for i in range(layout.ndim): - dl[:, i] = layout.dl[i] - - particles = Particles( - icells=h5_patch_grp["iCell"], - deltas=h5_patch_grp["delta"], - v=v, - weights=h5_patch_grp["weight"], - charges=h5_patch_grp["charge"], - dl=dl, - ) - - pdname = particle_dataset_name(basename) - if pdname in patch_datas: - raise ValueError("error - {} already in patchdata".format(pdname)) - - patch_datas[pdname] = ParticleData(layout, particles, pop_name(basename)) - - else: - for dataset_name in h5_patch_grp.keys(): - dataset = h5_patch_grp[dataset_name] - - if dataset_name not in field_qties: - raise RuntimeError( - "invalid dataset name : {} is not in {}".format( - dataset_name, field_qties - ) - ) - - pdata = FieldData(layout, field_qties[dataset_name], dataset) - - pdata_name = field_qties[dataset_name] - - if is_pop_fluid_file(basename): - pdata_name = pop_name(basename) + "_" + pdata_name - - if dataset_name in patch_datas: - raise ValueError("error - {} already in patchdata".format(dataset_name)) - - patch_datas[pdata_name] = pdata - - return True # valid patch assumed - - -def patch_has_datasets(h5_patch_grp): - return len(h5_patch_grp.keys()) > 0 - - -h5_time_grp_key = "t" - - -def hierarchy_fromh5(h5_filename, time, hier, silent=True): - import h5py - - data_file = h5py.File(h5_filename, "r") - basename = os.path.basename(h5_filename) - - root_cell_width = np.asarray(data_file.attrs["cell_width"]) - interp = data_file.attrs["interpOrder"] - domain_box = Box( - [0] * len(data_file.attrs["domain_box"]), data_file.attrs["domain_box"] - ) - - if create_from_all_times(time, hier): - # first create from first time - # then add all other times - if not silent: - print("creating hierarchy from all times in file") - times = list(data_file[h5_time_grp_key].keys()) - hier = hierarchy_fromh5(h5_filename, time=times[0], hier=hier, silent=silent) - if len(times) > 1: - for t in times[1:]: - hierarchy_fromh5(h5_filename, t, hier, silent=silent) - return hier - - if create_from_one_time(time, hier): - if not silent: - print("creating hierarchy from time {}".format(time)) - t = time - - h5_time_grp = data_file[h5_time_grp_key][time] - patch_levels = {} - - for plvl_key in h5_time_grp.keys(): - h5_patch_lvl_grp = h5_time_grp[plvl_key] - ilvl = int(plvl_key[2:]) - lvl_cell_width = root_cell_width / refinement_ratio**ilvl - patches = {} - - if ilvl not in patches: - patches[ilvl] = [] - - for pkey in h5_patch_lvl_grp.keys(): - h5_patch_grp = h5_patch_lvl_grp[pkey] - - patch_datas = {} - layout = make_layout(h5_patch_grp, lvl_cell_width, interp) - if patch_has_datasets(h5_patch_grp): - add_to_patchdata(patch_datas, h5_patch_grp, basename, layout) - - patches[ilvl].append( - Patch( - patch_datas, - h5_patch_grp.name.split("/")[-1], - layout=layout, - attrs={k: v for k, v in h5_patch_grp.attrs.items()}, - ) - ) - - patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl]) - - diag_hier = PatchHierarchy( - patch_levels, domain_box, refinement_ratio, t, data_file - ) - - return diag_hier - - if load_one_time(time, hier): - if not silent: - print("loading data at time {} into existing hierarchy".format(time)) - h5_time_grp = data_file[h5_time_grp_key][time] - t = time - - if t in hier.time_hier: - if not silent: - print("time already exist, adding data...") - - # time already exists in the hierarchy - # all we need to do is adding the data - # as patchDatas in the appropriate patches - # and levels, if data compatible with hierarchy - - patch_levels = hier.time_hier[t] - - for plvl_key in h5_time_grp.keys(): - ilvl = int(plvl_key[2:]) - lvl_cell_width = root_cell_width / refinement_ratio**ilvl - - for ipatch, pkey in enumerate(h5_time_grp[plvl_key].keys()): - h5_patch_grp = h5_time_grp[plvl_key][pkey] - - if patch_has_datasets(h5_patch_grp): - hier_patch = patch_levels[ilvl].patches[ipatch] - origin = h5_time_grp[plvl_key][pkey].attrs["origin"] - upper = h5_time_grp[plvl_key][pkey].attrs["upper"] - lower = h5_time_grp[plvl_key][pkey].attrs["lower"] - file_patch_box = Box(lower, upper) - - assert file_patch_box == hier_patch.box - assert (abs(origin - hier_patch.origin) < 1e-6).all() - assert (abs(lvl_cell_width - hier_patch.dl) < 1e-6).all() - - layout = make_layout(h5_patch_grp, lvl_cell_width, interp) - add_to_patchdata( - hier_patch.patch_datas, h5_patch_grp, basename, layout - ) - - return hier - - if not silent: - print("adding data to new time") - # time does not exist in the hierarchy - # we have to create a brand new set of patchLevels - # containing patches, and load data in their patchdatas - - patch_levels = {} - - for plvl_key in h5_time_grp.keys(): - ilvl = int(plvl_key[2:]) - - lvl_cell_width = root_cell_width / refinement_ratio**ilvl - lvl_patches = [] - - for ipatch, pkey in enumerate(h5_time_grp[plvl_key].keys()): - h5_patch_grp = h5_time_grp[plvl_key][pkey] - - layout = make_layout(h5_patch_grp, lvl_cell_width, interp) - patch_datas = {} - if patch_has_datasets(h5_patch_grp): - add_to_patchdata(patch_datas, h5_patch_grp, basename, layout) - lvl_patches.append( - Patch( - patch_datas, - h5_patch_grp.name.split("/")[-1], - layout=layout, - attrs={k: v for k, v in h5_patch_grp.attrs.items()}, - ) - ) - - patch_levels[ilvl] = PatchLevel(ilvl, lvl_patches) - - hier.time_hier[t] = patch_levels - return hier - - if load_all_times(time, hier): - if not silent: - print("loading all times in existing hier") - for time in data_file[h5_time_grp_key].keys(): - hier = hierarchy_fromh5(h5_filename, time, hier, silent=silent) - - return hier - - -def quantidic(ilvl, wrangler): - pl = wrangler.getPatchLevel(ilvl) - - return { - "density": pl.getDensity, - "bulkVelocity_x": pl.getVix, - "bulkVelocity_y": pl.getViy, - "bulkVelocity_z": pl.getViz, - "EM_B_x": pl.getBx, - "EM_B_y": pl.getBy, - "EM_B_z": pl.getBz, - "EM_E_x": pl.getEx, - "EM_E_y": pl.getEy, - "EM_E_z": pl.getEz, - "flux_x": pl.getFx, - "flux_y": pl.getFy, - "flux_z": pl.getFz, - "particles": pl.getParticles, - } - - -def isFieldQty(qty): - return qty in ( - "density", - "bulkVelocity_x", - "bulkVelocity_y", - "bulkVelocity_z", - "EM_B_x", - "EM_B_y", - "EM_B_z", - "EM_E_x", - "EM_E_y", - "EM_E_z", - "flux_x", - "flux_y", - "flux_z", - "momentumTensor_xx", - "momentumTensor_xy", - "momentumTensor_xz", - "momentumTensor_yy", - "momentumTensor_yz", - "momentumTensor_zz", - ) - - -def hierarchy_from_sim(simulator, qty, pop=""): - dw = simulator.data_wrangler() - nbr_levels = dw.getNumberOfLevels() - patch_levels = {} - - root_cell_width = simulator.cell_width() - domain_box = Box([0] * len(root_cell_width), simulator.domain_box()) - assert len(domain_box.ndim) == len(simulator.domain_box().ndim) - - for ilvl in range(nbr_levels): - lvl_cell_width = root_cell_width / refinement_ratio**ilvl - - patches = {ilvl: [] for ilvl in range(nbr_levels)} - getters = quantidic(ilvl, dw) - - if isFieldQty(qty): - wpatches = getters[qty]() - for patch in wpatches: - patch_datas = {} - lower = patch.lower - upper = patch.upper - origin = patch.origin - layout = GridLayout( - Box(lower, upper), - origin, - lvl_cell_width, - interp_order=simulator.interporder(), - ) - pdata = FieldData(layout, field_qties[qty], patch.data) - patch_datas[qty] = pdata - patches[ilvl].append(Patch(patch_datas)) - - elif qty == "particles": - if pop == "": - raise ValueError("must specify pop argument for particles") - # here the getter returns a dict like this - # {'protons': {'patchGhost': [, - # ], - # 'domain': [, - # ]}} - - # domain particles are assumed to always be here - # but patchGhost and levelGhost may not be, depending on the level - - populationdict = getters[qty](pop)[pop] - - dom_dw_patches = populationdict["domain"] - for patch in dom_dw_patches: - patch_datas = {} - - lower = patch.lower - upper = patch.upper - origin = patch.origin - layout = GridLayout( - Box(lower, upper), - origin, - lvl_cell_width, - interp_order=simulator.interp_order(), - ) - v = np.asarray(patch.data.v).reshape(int(len(patch.data.v) / 3), 3) - - domain_particles = Particles( - icells=np.asarray(patch.data.iCell), - deltas=np.asarray(patch.data.delta), - v=v, - weights=np.asarray(patch.data.weight), - charges=np.asarray(patch.data.charge), - ) - - patch_datas[pop + "_particles"] = ParticleData( - layout, domain_particles, pop - ) - patches[ilvl].append(Patch(patch_datas)) - - # ok now let's add the patchGhost if present - # note that patchGhost patches may not be the same list as the - # domain patches... since not all patches may not have patchGhost while they do have - # domain... while looping on the patchGhost items, we need to search in - # the already created patches which one to which add the patchGhost particles - - for ghostParticles in ["patchGhost", "levelGhost"]: - if ghostParticles in populationdict: - for dwpatch in populationdict[ghostParticles]: - v = np.asarray(dwpatch.data.v) - s = v.size - v = v[:].reshape(int(s / 3), 3) - - patchGhost_part = Particles( - icells=np.asarray(dwpatch.data.iCell), - deltas=np.asarray(dwpatch.data.delta), - v=v, - weights=np.asarray(dwpatch.data.weight), - charges=np.asarray(dwpatch.data.charge), - ) - - box = Box(dwpatch.lower, dwpatch.upper) - - # now search which of the already created patches has the same box - # once found we add the new particles to the ones already present - - patch = [p for p in patches[ilvl] if p.box == box][0] - patch.patch_datas[pop + "_particles"].dataset.add( - patchGhost_part - ) - - else: - raise ValueError("{} is not a valid quantity".format(qty)) - - patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl]) - - return PatchHierarchy(patch_levels, domain_box, time=simulator.currentTime()) - - -def hierarchy_from( - simulator=None, qty=None, pop="", h5_filename=None, time=None, hier=None -): - """ - this function reads an HDF5 PHARE file and returns a PatchHierarchy from - which data is accessible. - if 'time' is None, all times in the file will be read, if a time is given - then only that time will be read - if 'hier' is None, then a new hierarchy will be created, if not then the - given hierarchy 'hier' will be filled. - - The function fails if the data is already in hierarchy - """ - - if simulator is not None and h5_filename is not None: - raise ValueError("cannot pass both a simulator and a h5 file") - - if h5_filename is not None: - return hierarchy_fromh5(h5_filename, time, hier) - - if simulator is not None and qty is not None: - return hierarchy_from_sim(simulator, qty, pop=pop) - - raise ValueError("can't make hierarchy") - - -def merge_particles(hierarchy): - for time, patch_levels in hierarchy.time_hier.items(): - for ilvl, plvl in patch_levels.items(): - for ip, patch in enumerate(plvl.patches): - pdatas = patch.patch_datas - domain_pdata = [ - (pdname, pd) for pdname, pd in pdatas.items() if "domain" in pdname - ][0] - - pghost_pdatas = [ - (pdname, pd) - for pdname, pd in pdatas.items() - if "patchGhost" in pdname - ] - lghost_pdatas = [ - (pdname, pd) - for pdname, pd in pdatas.items() - if "levelGhost" in pdname - ] - - pghost_pdata = pghost_pdatas[0] if pghost_pdatas else None - lghost_pdata = lghost_pdatas[0] if lghost_pdatas else None - - if pghost_pdata is not None: - domain_pdata[1].dataset.add(pghost_pdata[1].dataset) - del pdatas[pghost_pdata[0]] - - if lghost_pdata is not None: - domain_pdata[1].dataset.add(lghost_pdata[1].dataset) - del pdatas[lghost_pdata[0]] - - popname = domain_pdata[0].split("_")[0] - pdatas[popname + "_particles"] = pdatas[domain_pdata[0]] - del pdatas[domain_pdata[0]] - - -def h5_filename_from(diagInfo): - # diagInfo.quantity starts with a / , hence [1:] - return (diagInfo.quantity + ".h5").replace("/", "_")[1:] - - -def get_times_from_h5(filepath): - import h5py - - f = h5py.File(filepath, "r") - times = np.array(sorted([float(s) for s in list(f["t"].keys())])) - f.close() - return times - - -def getPatch(hier, point): - """ - convenience function mainly for debug. returns a dict - {ilevel:patch} for patches in which the given point is - """ - patches = {} - counts = {ilvl: 0 for ilvl in hier.levels().keys()} - for ilvl, lvl in hier.levels().items(): - for p in lvl.patches: - px, py = point - dx, dy = p.layout.dl - ix = int(px / dx) - iy = int(py / dy) - if (ix, iy) in p.box: - patches[ilvl] = p - counts[ilvl] += 1 - for k, v in counts.items(): - if v > 1: - print("error : ", k, v) - raise RuntimeError("more than one patch found for point") - return patches diff --git a/pyphare/pyphare/pharesee/hierarchy/__init__.py b/pyphare/pyphare/pharesee/hierarchy/__init__.py new file mode 100644 index 000000000..edeed0416 --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/__init__.py @@ -0,0 +1,32 @@ +from .hierarchy import PatchHierarchy +from .fromh5 import hierarchy_fromh5 +from .fromsim import hierarchy_from_sim + + +def hierarchy_from( + simulator=None, qty=None, pop="", h5_filename=None, time=None, hier=None +): + from .fromh5 import hierarchy_fromh5 + from .fromsim import hierarchy_from_sim + + """ + this function reads an HDF5 PHARE file and returns a PatchHierarchy from + which data is accessible. + if 'time' is None, all times in the file will be read, if a time is given + then only that time will be read + if 'hier' is None, then a new hierarchy will be created, if not then the + given hierarchy 'hier' will be filled. + + The function fails if the data is already in hierarchy + """ + + if simulator is not None and h5_filename is not None: + raise ValueError("cannot pass both a simulator and a h5 file") + + if h5_filename is not None: + return hierarchy_fromh5(h5_filename, time, hier) + + if simulator is not None and qty is not None: + return hierarchy_from_sim(simulator, qty, pop=pop) + + raise ValueError("can't make hierarchy") diff --git a/pyphare/pyphare/pharesee/hierarchy/fromh5.py b/pyphare/pyphare/pharesee/hierarchy/fromh5.py new file mode 100644 index 000000000..710a82c2e --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/fromh5.py @@ -0,0 +1,306 @@ +import os +import numpy as np + +from .patch import Patch +from .patchlevel import PatchLevel +from .patchdata import FieldData, ParticleData +from ..particles import Particles +from .hierarchy import PatchHierarchy +from ...core.box import Box +from ...core.phare_utilities import refinement_ratio +from ...core.gridlayout import GridLayout +from .hierarchy_utils import field_qties + + +particle_files_patterns = ("domain", "patchGhost", "levelGhost") + + +def make_layout(h5_patch_grp, cell_width, interp_order): + origin = h5_patch_grp.attrs["origin"] + upper = h5_patch_grp.attrs["upper"] + lower = h5_patch_grp.attrs["lower"] + return GridLayout(Box(lower, upper), origin, cell_width, interp_order=interp_order) + + +def is_pop_fluid_file(basename): + return (is_particle_file(basename) is False) and "pop" in basename + + +def particle_dataset_name(basename): + """ + return "alpha_domain" from ions_pop_alpha_domain.h5 + """ + popname = basename.strip(".h5").split("_")[-2] + part_type = basename.strip(".h5").split("_")[-1] + dataset_name = popname + "_" + part_type + + return dataset_name + + +def is_particle_file(filename): + return any([pattern in filename for pattern in particle_files_patterns]) + + +def pop_name(basename): + return basename.strip(".h5").split("_")[2] + + +def add_to_patchdata(patch_datas, h5_patch_grp, basename, layout): + """ + adds data in the h5_patch_grp in the given PatchData dict + returns True if valid h5 patch found + """ + + if is_particle_file(basename): + v = np.asarray(h5_patch_grp["v"]) + s = v.size + v = v[:].reshape(int(s / 3), 3) + nbrParts = v.shape[0] + dl = np.zeros((nbrParts, layout.ndim)) + for i in range(layout.ndim): + dl[:, i] = layout.dl[i] + + particles = Particles( + icells=h5_patch_grp["iCell"], + deltas=h5_patch_grp["delta"], + v=v, + weights=h5_patch_grp["weight"], + charges=h5_patch_grp["charge"], + dl=dl, + ) + + pdname = particle_dataset_name(basename) + if pdname in patch_datas: + raise ValueError("error - {} already in patchdata".format(pdname)) + + patch_datas[pdname] = ParticleData(layout, particles, pop_name(basename)) + + else: + for dataset_name in h5_patch_grp.keys(): + dataset = h5_patch_grp[dataset_name] + + if dataset_name not in field_qties: + raise RuntimeError( + "invalid dataset name : {} is not in {}".format( + dataset_name, field_qties + ) + ) + + pdata = FieldData(layout, field_qties[dataset_name], dataset) + + pdata_name = field_qties[dataset_name] + + if is_pop_fluid_file(basename): + pdata_name = pop_name(basename) + "_" + pdata_name + + if dataset_name in patch_datas: + raise ValueError("error - {} already in patchdata".format(dataset_name)) + + patch_datas[pdata_name] = pdata + + return True # valid patch assumed + + +def patch_has_datasets(h5_patch_grp): + return len(h5_patch_grp.keys()) > 0 + + +def h5_filename_from(diagInfo): + # diagInfo.quantity starts with a / , hence [1:] + return (diagInfo.quantity + ".h5").replace("/", "_")[1:] + + +def get_times_from_h5(filepath): + import h5py + + f = h5py.File(filepath, "r") + times = np.array(sorted([float(s) for s in list(f["t"].keys())])) + f.close() + return times + + +h5_time_grp_key = "t" + + +def create_from_all_times(time, hier): + return time is None and hier is None + + +def create_from_one_time(time, hier): + return time is not None and hier is None + + +def load_all_times(time, hier): + return time is None and hier is not None + + +def load_one_time(time, hier): + return time is not None and hier is not None + + +def fromh5_single_time(h5_filename, time, silent=True): + pass + + +def hierarchy_fromh5_(h5_filename, time, hier, silent=True): + + import h5py + + data_file = h5py.File(h5_filename, "r") + basename = os.path.basename(h5_filename) + root_cell_width = np.asarray(data_file.attrs["cell_width"]) + interp = data_file.attrs["interpOrder"] + dimension = len(data_file.attrs["domain_box"]) + domain_box = Box([0] * dimension, data_file.attrs["domain_box"]) + + +def hierarchy_fromh5(h5_filename, time, hier, silent=True): + import h5py + + data_file = h5py.File(h5_filename, "r") + basename = os.path.basename(h5_filename) + + root_cell_width = np.asarray(data_file.attrs["cell_width"]) + interp = data_file.attrs["interpOrder"] + domain_box = Box( + [0] * len(data_file.attrs["domain_box"]), data_file.attrs["domain_box"] + ) + + if create_from_all_times(time, hier): + # first create from first time + # then add all other times + if not silent: + print("creating hierarchy from all times in file") + times = list(data_file[h5_time_grp_key].keys()) + hier = hierarchy_fromh5(h5_filename, time=times[0], hier=hier, silent=silent) + if len(times) > 1: + for t in times[1:]: + hierarchy_fromh5(h5_filename, t, hier, silent=silent) + return hier + + if create_from_one_time(time, hier): + if not silent: + print("creating hierarchy from time {}".format(time)) + t = time + + h5_time_grp = data_file[h5_time_grp_key][time] + patch_levels = {} + + for plvl_key in h5_time_grp.keys(): + h5_patch_lvl_grp = h5_time_grp[plvl_key] + ilvl = int(plvl_key[2:]) + lvl_cell_width = root_cell_width / refinement_ratio**ilvl + patches = {} + + if ilvl not in patches: + patches[ilvl] = [] + + for pkey in h5_patch_lvl_grp.keys(): + h5_patch_grp = h5_patch_lvl_grp[pkey] + + patch_datas = {} + layout = make_layout(h5_patch_grp, lvl_cell_width, interp) + if patch_has_datasets(h5_patch_grp): + add_to_patchdata(patch_datas, h5_patch_grp, basename, layout) + + patches[ilvl].append( + Patch( + patch_datas, + h5_patch_grp.name.split("/")[-1], + layout=layout, + attrs={k: v for k, v in h5_patch_grp.attrs.items()}, + ) + ) + + patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl]) + + diag_hier = PatchHierarchy( + patch_levels, domain_box, refinement_ratio, t, data_file + ) + + return diag_hier + + if load_one_time(time, hier): + if not silent: + print("loading data at time {} into existing hierarchy".format(time)) + h5_time_grp = data_file[h5_time_grp_key][time] + t = time + + if t in hier.time_hier: + if not silent: + print("time already exist, adding data...") + + # time already exists in the hierarchy + # all we need to do is adding the data + # as patchDatas in the appropriate patches + # and levels, if data compatible with hierarchy + + patch_levels = hier.time_hier[t] + + for plvl_key in h5_time_grp.keys(): + ilvl = int(plvl_key[2:]) + lvl_cell_width = root_cell_width / refinement_ratio**ilvl + + for ipatch, pkey in enumerate(h5_time_grp[plvl_key].keys()): + h5_patch_grp = h5_time_grp[plvl_key][pkey] + + if patch_has_datasets(h5_patch_grp): + hier_patch = patch_levels[ilvl].patches[ipatch] + origin = h5_time_grp[plvl_key][pkey].attrs["origin"] + upper = h5_time_grp[plvl_key][pkey].attrs["upper"] + lower = h5_time_grp[plvl_key][pkey].attrs["lower"] + file_patch_box = Box(lower, upper) + + assert file_patch_box == hier_patch.box + assert (abs(origin - hier_patch.origin) < 1e-6).all() + assert (abs(lvl_cell_width - hier_patch.dl) < 1e-6).all() + + layout = make_layout(h5_patch_grp, lvl_cell_width, interp) + add_to_patchdata( + hier_patch.patch_datas, h5_patch_grp, basename, layout + ) + + return hier + + if not silent: + print("adding data to new time") + # time does not exist in the hierarchy + # we have to create a brand new set of patchLevels + # containing patches, and load data in their patchdatas + + patch_levels = {} + + for plvl_key in h5_time_grp.keys(): + ilvl = int(plvl_key[2:]) + + lvl_cell_width = root_cell_width / refinement_ratio**ilvl + lvl_patches = [] + + for ipatch, pkey in enumerate(h5_time_grp[plvl_key].keys()): + h5_patch_grp = h5_time_grp[plvl_key][pkey] + + layout = make_layout(h5_patch_grp, lvl_cell_width, interp) + patch_datas = {} + if patch_has_datasets(h5_patch_grp): + add_to_patchdata(patch_datas, h5_patch_grp, basename, layout) + lvl_patches.append( + Patch( + patch_datas, + h5_patch_grp.name.split("/")[-1], + layout=layout, + attrs={k: v for k, v in h5_patch_grp.attrs.items()}, + ) + ) + + patch_levels[ilvl] = PatchLevel(ilvl, lvl_patches) + + hier.time_hier[t] = patch_levels + return hier + + if load_all_times(time, hier): + if not silent: + print("loading all times in existing hier") + for time in data_file[h5_time_grp_key].keys(): + hier = hierarchy_fromh5(h5_filename, time, hier, silent=silent) + + return hier diff --git a/pyphare/pyphare/pharesee/hierarchy/fromsim.py b/pyphare/pyphare/pharesee/hierarchy/fromsim.py new file mode 100644 index 000000000..ecdc24600 --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/fromsim.py @@ -0,0 +1,123 @@ +from .hierarchy_utils import isFieldQty, field_qties, quantidic, refinement_ratio +from .patchdata import FieldData, ParticleData +from .patch import Patch +from .patchlevel import PatchLevel +from .hierarchy import PatchHierarchy +from ..particles import Particles +from ...core.gridlayout import GridLayout +from ...core.box import Box + +import numpy as np + + +def hierarchy_from_sim(simulator, qty, pop=""): + dw = simulator.data_wrangler() + nbr_levels = dw.getNumberOfLevels() + patch_levels = {} + + root_cell_width = simulator.cell_width() + domain_box = Box([0] * len(root_cell_width), simulator.domain_box()) + assert len(domain_box.ndim) == len(simulator.domain_box().ndim) + + for ilvl in range(nbr_levels): + lvl_cell_width = root_cell_width / refinement_ratio**ilvl + + patches = {ilvl: [] for ilvl in range(nbr_levels)} + getters = quantidic(ilvl, dw) + + if isFieldQty(qty): + wpatches = getters[qty]() + for patch in wpatches: + patch_datas = {} + lower = patch.lower + upper = patch.upper + origin = patch.origin + layout = GridLayout( + Box(lower, upper), + origin, + lvl_cell_width, + interp_order=simulator.interporder(), + ) + pdata = FieldData(layout, field_qties[qty], patch.data) + patch_datas[qty] = pdata + patches[ilvl].append(Patch(patch_datas)) + + elif qty == "particles": + if pop == "": + raise ValueError("must specify pop argument for particles") + # here the getter returns a dict like this + # {'protons': {'patchGhost': [, + # ], + # 'domain': [, + # ]}} + + # domain particles are assumed to always be here + # but patchGhost and levelGhost may not be, depending on the level + + populationdict = getters[qty](pop)[pop] + + dom_dw_patches = populationdict["domain"] + for patch in dom_dw_patches: + patch_datas = {} + + lower = patch.lower + upper = patch.upper + origin = patch.origin + layout = GridLayout( + Box(lower, upper), + origin, + lvl_cell_width, + interp_order=simulator.interp_order(), + ) + v = np.asarray(patch.data.v).reshape(int(len(patch.data.v) / 3), 3) + + domain_particles = Particles( + icells=np.asarray(patch.data.iCell), + deltas=np.asarray(patch.data.delta), + v=v, + weights=np.asarray(patch.data.weight), + charges=np.asarray(patch.data.charge), + ) + + patch_datas[pop + "_particles"] = ParticleData( + layout, domain_particles, pop + ) + patches[ilvl].append(Patch(patch_datas)) + + # ok now let's add the patchGhost if present + # note that patchGhost patches may not be the same list as the + # domain patches... since not all patches may not have patchGhost while they do have + # domain... while looping on the patchGhost items, we need to search in + # the already created patches which one to which add the patchGhost particles + + for ghostParticles in ["patchGhost", "levelGhost"]: + if ghostParticles in populationdict: + for dwpatch in populationdict[ghostParticles]: + v = np.asarray(dwpatch.data.v) + s = v.size + v = v[:].reshape(int(s / 3), 3) + + patchGhost_part = Particles( + icells=np.asarray(dwpatch.data.iCell), + deltas=np.asarray(dwpatch.data.delta), + v=v, + weights=np.asarray(dwpatch.data.weight), + charges=np.asarray(dwpatch.data.charge), + ) + + box = Box(dwpatch.lower, dwpatch.upper) + + # now search which of the already created patches has the same box + # once found we add the new particles to the ones already present + + patch = [p for p in patches[ilvl] if p.box == box][0] + patch.patch_datas[pop + "_particles"].dataset.add( + patchGhost_part + ) + + else: + raise ValueError("{} is not a valid quantity".format(qty)) + + patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl]) + + return PatchHierarchy(patch_levels, domain_box, time=simulator.currentTime()) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py new file mode 100644 index 000000000..a1713dd3b --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -0,0 +1,674 @@ +from .patch import Patch +from .patchlevel import PatchLevel +from ...core.box import Box +from ...core import box as boxm +from ...core.phare_utilities import refinement_ratio + +import numpy as np +import matplotlib.pyplot as plt + + +class PatchHierarchy(object): + """is a collection of patch levels""" + + def __init__( + self, + patch_levels, + domain_box, + refinement_ratio=2, + time=0.0, + data_files=None, + **kwargs, + ): + self.patch_levels = patch_levels + self.ndim = len(domain_box.lower) + self.time_hier = {} + self.time_hier.update({self.format_timestamp(time): patch_levels}) + + self.domain_box = domain_box + self.refinement_ratio = refinement_ratio + + self.data_files = {} + self._sim = None + + if data_files is not None: + self.data_files.update(data_files) + + if len(self.quantities()) > 1: + for qty in self.quantities(): + if qty in self.__dict__: + continue + first = True + for time, levels in self.time_hier.items(): + new_lvls = {} + for ilvl, level in levels.items(): + patches = [] + for patch in level.patches: + patches += [Patch({qty: patch.patch_datas[qty]}, patch.id)] + new_lvls[ilvl] = PatchLevel(ilvl, patches) + if first: + self.__dict__[qty] = PatchHierarchy( + new_lvls, self.domain_box, time=time + ) + first = False + else: + self.qty.time_hier[time] = new_lvls # pylint: disable=E1101 + + def __getitem__(self, qty): + return self.__dict__[qty] + + @property + def sim(self): + if self._sim: + return self._sim + + if "py_attrs" not in self.data_files: + raise ValueError("Simulation is not available for deserialization") + + from ...pharein.simulation import deserialize + + try: + self._sim = deserialize( + self.data_files["py_attrs"].attrs["serialized_simulation"] + ) + except Exception as e: + raise RuntimeError(f"Failed to deserialize simulation from data file : {e}") + return self._sim + + def __call__(self, qty=None, **kwargs): + # take slice/slab of 1/2d array from 2/3d array + def cuts(c, coord): + return c > coord.min() and c < coord.max() + + class Extractor: + def __init__(self): + self.exclusions = [] + + def extract(self, coord, data): + mask = coord == coord + for exclusion in self.exclusions: + idx = np.where( + (coord > exclusion[0] - 1e-6) & (coord < exclusion[1] + 1e-6) + )[0] + mask[idx] = False + + self.exclusions += [(coord.min(), coord.max())] + return coord[mask], data[mask] + + def domain_coords(patch, qty): + pd = patch.patch_datas[qty] + nbrGhosts = pd.ghosts_nbr[0] + return pd.x[nbrGhosts:-nbrGhosts], pd.y[nbrGhosts:-nbrGhosts] + + if len(kwargs) < 1 or len(kwargs) > 3: + raise ValueError("Error - must provide coordinates") + if qty is None: + if len(self.quantities()) == 1: + qty = self.quantities()[0] + else: + raise ValueError( + "The PatchHierarchy has several quantities but none is specified" + ) + + if len(kwargs) == 1: # 1D cut + if "x" in kwargs: + c = kwargs["x"] + slice_dim = 1 + cst_dim = 0 + else: + c = kwargs["y"] + slice_dim = 0 + cst_dim = 1 + + extractor = Extractor() + datas = [] + coords = [] + ilvls = list(self.levels().keys())[::-1] + + for ilvl in ilvls: + lvl = self.patch_levels[ilvl] + for patch in lvl.patches: + slice_coord = domain_coords(patch, qty)[slice_dim] + cst_coord = domain_coords(patch, qty)[cst_dim] + + if cuts(c, cst_coord): + data = patch(qty, **kwargs) + coord_keep, data_keep = extractor.extract(slice_coord, data) + datas += [data_keep] + coords += [coord_keep] + + cut = np.concatenate(datas) + coords = np.concatenate(coords) + ic = np.argsort(coords) + coords = coords[ic] + cut = cut[ic] + return coords, cut + + def _default_time(self): + return self.times()[0] + + def finest_level(self, time=None): + if time is None: + time = self._default_time() + return max(list(self.levels(time=time).keys())) + + def levels(self, time=None): + if time is None: + time = self._default_time() + return self.time_hier[self.format_timestamp(time)] + + def level(self, level_number, time=None): + return self.levels(time)[level_number] + + def levelNbr(self, time=None): + if time is None: + time = self._default_time() + return len(self.levels(time).items()) + + def levelNbrs(self, time=None): + if time is None: + time = self._default_time() + return list(self.levels(time).keys()) + + def is_homogeneous(self): + """ + return True if all patches of all levels at all times + have the same patch data quantities + """ + qties = self._quantities() + it_is = True + for time, levels in self.time_hier.items(): + for ilvl, lvl in levels.items(): + for patch in lvl.patches: + pdnames = list(patch.patch_datas.keys()) + if len(pdnames): # do not compare empty patches + it_is &= qties == pdnames + return it_is + + def _quantities(self): + # we return the name of the patchdatas of the first level that has + # patches with data. checking that patchdatas are not {} is important + # since some patches might be empty (e.g. level ghost patchdatas on L0) + for ilvl, lvl in self.levels().items(): + if len(lvl.patches) > 0 and len(lvl.patches[0].patch_datas): + return list(self.level(ilvl).patches[0].patch_datas.keys()) + return [] + + def quantities(self): + if not self.is_homogeneous(): + raise RuntimeError("Error - hierarchy is not homogeneous") + return self._quantities() + + def global_min(self, qty, **kwargs): + time = kwargs.get("time", self._default_time()) + first = True + for ilvl, lvl in self.levels(time).items(): + for patch in lvl.patches: + pd = patch.patch_datas[qty] + if first: + m = pd.dataset[:].min() + first = False + else: + m = min(m, pd.dataset[:].min()) + + return m + + def global_max(self, qty, **kwargs): + time = kwargs.get("time", self._default_time()) + first = True + for _, lvl in self.levels(time).items(): + for patch in lvl.patches: + pd = patch.patch_datas[qty] + if first: + m = pd.dataset[:].max() + first = False + else: + m = max(m, pd.dataset[:].max()) + + return m + + def refined_domain_box(self, level_number): + """ + returns the domain box refined for a given level number + """ + assert level_number >= 0 + return boxm.refine(self.domain_box, self.refinement_ratio**level_number) + + def format_timestamp(self, timestamp): + if isinstance(timestamp, str): + return timestamp + return "{:.10f}".format(timestamp) + + def level_domain_box(self, level_number): + if level_number == 0: + return self.domain_box + return self.refined_domain_box(level_number) + + def __str__(self): + s = "Hierarchy: \n" + for t, patch_levels in self.time_hier.items(): + for ilvl, lvl in patch_levels.items(): + s = s + "Level {}\n".format(ilvl) + for ip, patch in enumerate(lvl.patches): + for qty_name, pd in patch.patch_datas.items(): + pdstr = " P{ip} {type} {pdname} box is {box} and ghost box is {gbox}" + s = s + pdstr.format( + ip=ip, + type=type(pd.dataset), + pdname=qty_name, + box=patch.box, + gbox=pd.ghost_box, + ) + s = s + "\n" + return s + + def times(self): + return np.sort(np.asarray(list(self.time_hier.keys()))) + + def plot_patches(self, save=False): + fig, ax = plt.subplots(figsize=(10, 3)) + for ilvl, lvl in self.levels(0.0).items(): + lvl_offset = ilvl * 0.1 + for patch in lvl.patches: + dx = patch.dl[0] + x0 = patch.box.lower * dx + x1 = patch.box.upper * dx + xcells = np.arange(x0, x1 + dx, dx) + y = lvl_offset + np.zeros_like(xcells) + ax.plot(xcells, y, marker=".") + + if save: + fig.savefig("hierarchy.png") + + def box_to_Rectangle(self, box): + from matplotlib.patches import Rectangle + + return Rectangle(box.lower, *box.shape) + + def plot_2d_patches(self, ilvl, collections, **kwargs): + if isinstance(collections, list) and all( + [isinstance(el, Box) for el in collections] + ): + collections = [{"boxes": collections}] + + from matplotlib.collections import PatchCollection + + level_domain_box = self.level_domain_box(ilvl) + mi, ma = level_domain_box.lower.min(), level_domain_box.upper.max() + + fig, ax = kwargs.get("subplot", plt.subplots(figsize=(6, 6))) + + for collection in collections: + facecolor = collection.get("facecolor", "none") + edgecolor = collection.get("edgecolor", "purple") + alpha = collection.get("alpha", 1) + rects = [self.box_to_Rectangle(box) for box in collection["boxes"]] + + ax.add_collection( + PatchCollection( + rects, facecolor=facecolor, alpha=alpha, edgecolor=edgecolor + ) + ) + + if "title" in kwargs: + from textwrap import wrap + + xfigsize = int(fig.get_size_inches()[0] * 10) # 10 characters per inch + ax.set_title("\n".join(wrap(kwargs["title"], xfigsize))) + + major_ticks = np.arange(mi - 5, ma + 5 + 5, 5) + ax.set_xticks(major_ticks) + ax.set_yticks(major_ticks) + + minor_ticks = np.arange(mi - 5, ma + 5 + 5, 1) + ax.set_xticks(minor_ticks, minor=True) + ax.set_yticks(minor_ticks, minor=True) + + ax.grid(which="both") + + return fig + + def plot1d(self, **kwargs): + """ + plot + """ + usr_lvls = kwargs.get("levels", (0,)) + qty = kwargs.get("qty", None) + time = kwargs.get("time", self.times()[0]) + + if "ax" not in kwargs: + fig, ax = plt.subplots() + else: + ax = kwargs["ax"] + fig = ax.figure + for lvl_nbr, level in self.levels(time).items(): + if lvl_nbr not in usr_lvls: + continue + for ip, patch in enumerate(level.patches): + pdata_nbr = len(patch.patch_datas) + pdata_names = list(patch.patch_datas.keys()) + if qty is None and pdata_nbr != 1: + multiple = "multiple quantities in patch, " + err = ( + multiple + + "please specify a quantity in " + + " ".join(pdata_names) + ) + raise ValueError(err) + if qty is None: + qty = pdata_names[0] + + layout = patch.patch_datas[qty].layout + nbrGhosts = layout.nbrGhostFor(qty) + val = patch.patch_datas[qty][patch.box] + x = patch.patch_datas[qty].x[nbrGhosts[0] : -nbrGhosts[0]] + label = "L{level}P{patch}".format(level=lvl_nbr, patch=ip) + marker = kwargs.get("marker", "") + ls = kwargs.get("ls", "--") + color = kwargs.get("color", "k") + ax.plot(x, val, label=label, marker=marker, ls=ls, color=color) + + ax.set_title(kwargs.get("title", "")) + ax.set_xlabel(kwargs.get("xlabel", "x")) + ax.set_ylabel(kwargs.get("ylabel", qty)) + if "xlim" in kwargs: + ax.set_xlim(kwargs["xlim"]) + if "ylim" in kwargs: + ax.set_ylim(kwargs["ylim"]) + + if kwargs.get("legend", None) is not None: + ax.legend() + + if "filename" in kwargs: + fig.savefig(kwargs["filename"]) + + def plot2d(self, **kwargs): + from matplotlib.patches import Rectangle + from mpl_toolkits.axes_grid1 import make_axes_locatable + + time = kwargs.get("time", self._default_time()) + usr_lvls = kwargs.get("levels", self.levelNbrs(time)) + default_qty = None + if len(self.quantities()) == 1: + default_qty = self.quantities()[0] + qty = kwargs.get("qty", default_qty) + + if "ax" not in kwargs: + fig, ax = plt.subplots() + else: + ax = kwargs["ax"] + fig = ax.figure + + glob_min = self.global_min(qty) + glob_max = self.global_max(qty) + # assumes max 5 levels... + patchcolors = {ilvl: "k" for ilvl in usr_lvls} + patchcolors = kwargs.get("patchcolors", patchcolors) + if not isinstance(patchcolors, dict): + patchcolors = dict(zip(usr_lvls, patchcolors)) + + linewidths = {ilvl: 1 for ilvl in usr_lvls} + linewidths = kwargs.get("lw", linewidths) + if not isinstance(linewidths, dict): + linewidths = dict(zip(usr_lvls, linewidths)) + linestyles = {ilvl: "-" for ilvl in usr_lvls} + linestyles = kwargs.get("ls", linestyles) + if not isinstance(linestyles, dict): + linestyles = dict(zip(usr_lvls, linestyles)) + + for lvl_nbr, lvl in self.levels(time).items(): + if lvl_nbr not in usr_lvls: + continue + for patch in self.level(lvl_nbr, time).patches: + pdat = patch.patch_datas[qty] + data = pdat.dataset[:] + nbrGhosts = pdat.ghosts_nbr + x = pdat.x + y = pdat.y + + # if nbrGhosts is 0, we cannot do array[0,-0] + if np.all(nbrGhosts == np.zeros_like(nbrGhosts)): + x = np.copy(x) + y = np.copy(y) + else: + data = pdat[patch.box] + x = np.copy(x[nbrGhosts[0] : -nbrGhosts[0]]) + y = np.copy(y[nbrGhosts[1] : -nbrGhosts[1]]) + dx, dy = pdat.layout.dl + x -= dx * 0.5 + y -= dy * 0.5 + x = np.append(x, x[-1] + dx) + y = np.append(y, y[-1] + dy) + im = ax.pcolormesh( + x, + y, + data.T, + cmap=kwargs.get("cmap", "Spectral_r"), + vmin=kwargs.get("vmin", glob_min - 1e-6), + vmax=kwargs.get("vmax", glob_max + 1e-6), + ) + + if kwargs.get("plot_patches", False) is True: + r = Rectangle( + (patch.box.lower[0] * dx, patch.box.lower[1] * dy), + patch.box.shape[0] * dx, + patch.box.shape[1] * dy, + fc="none", + ec=patchcolors[lvl_nbr], + alpha=0.4, + lw=linewidths[lvl_nbr], + ls=linestyles[lvl_nbr], + ) + ax.add_patch(r) + + ax.set_aspect(kwargs.get("aspect", "equal")) + ax.set_title(kwargs.get("title", "")) + ax.set_xlabel(kwargs.get("xlabel", "x")) + ax.set_ylabel(kwargs.get("ylabel", "y")) + if "xlim" in kwargs: + ax.set_xlim(kwargs["xlim"]) + if "ylim" in kwargs: + ax.set_ylim(kwargs["ylim"]) + + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.08) + plt.colorbar(im, ax=ax, cax=cax) + + if kwargs.get("legend", None) is not None: + ax.legend() + + if "filename" in kwargs: + fig.savefig(kwargs["filename"], dpi=kwargs.get("dpi", 200)) + + return fig, ax + + def plot(self, **kwargs): + if self.ndim == 1: + return self.plot1d(**kwargs) + elif self.ndim == 2: + return self.plot2d(**kwargs) + + def dist_plot(self, **kwargs): + """ + plot phase space of a particle hierarchy + """ + import copy + + from ..plotting import dist_plot as dp + + usr_lvls = kwargs.get("levels", (0,)) + finest = kwargs.get("finest", False) + pops = kwargs.get("pop", []) + time = kwargs.get("time", self.times()[0]) + axis = kwargs.get("axis", ("Vx", "Vy")) + all_pops = list(self.level(0, time).patches[0].patch_datas.keys()) + + vmin = kwargs.get("vmin", -2) + vmax = kwargs.get("vmax", 2) + dv = kwargs.get("dv", 0.05) + vbins = vmin + dv * np.arange(int((vmax - vmin) / dv)) + + if finest: + final = finest_part_data(self) + if axis[0] == "x": + xbins = amr_grid(self, time) + bins = (xbins, vbins) + else: + bins = (vbins, vbins) + kwargs["bins"] = bins + + else: + final = {pop: None for pop in all_pops} + for lvl_nbr, level in self.levels(time).items(): + if lvl_nbr not in usr_lvls: + continue + for ip, patch in enumerate(level.patches): + if len(pops) == 0: + pops = list(patch.patch_datas.keys()) + + for pop in pops: + tmp = copy.copy(patch.patch_datas[pop].dataset) + + if final[pop] is None: + final[pop] = tmp + else: + final[pop].add(tmp) + + # select particles + if "select" in kwargs: + for pop, particles in final.items(): + final[pop] = kwargs["select"](particles) + + return final, dp(final, **kwargs) + + +def finest_part_data(hierarchy, time=None): + """ + returns a dict {popname : Particles} + Particles contained in the dict are those from + the finest patches available at a given location + """ + from copy import deepcopy + + from ..particles import remove + + # we are going to return a dict {popname : Particles} + # we prepare it with population names + aPatch = hierarchy.level(0, time=time).patches[0] + particles = {popname: None for popname in aPatch.patch_datas.keys()} + + # our strategy is to explore the hierarchy from the finest + # level to the coarsest. at Each level we keep only particles + # that are in cells that are not overlaped by finer boxes + + # this dict keeps boxes for patches at each level + # each level will thus need this dict to see next finer boxes + lvlPatchBoxes = {ilvl: [] for ilvl in range(hierarchy.finest_level(time) + 1)} + + for ilvl in range(hierarchy.finest_level(time) + 1)[::-1]: + plvl = hierarchy.level(ilvl, time=time) + for ip, patch in enumerate(plvl.patches): + lvlPatchBoxes[ilvl].append(patch.box) + for popname, pdata in patch.patch_datas.items(): + # if we're at the finest level + # we need to keep all particles + if ilvl == hierarchy.finest_level(time): + if particles[popname] is None: + particles[popname] = deepcopy(pdata.dataset) + else: + particles[popname].add(deepcopy(pdata.dataset)) + + # if there is a finer level + # we need to keep only those of the current patch + # that are not in cells overlaped by finer patch boxes + else: + icells = pdata.dataset.iCells + parts = deepcopy(pdata.dataset) + create = True + for finerBox in lvlPatchBoxes[ilvl + 1]: + coarseFinerBox = boxm.coarsen(finerBox, refinement_ratio) + within = np.where( + (icells >= coarseFinerBox.lower[0]) + & (icells <= coarseFinerBox.upper[0]) + )[0] + if create: + toRemove = within + create = False + else: + toRemove = np.concatenate((toRemove, within)) + + if toRemove.size != 0: + parts = remove(parts, toRemove) + if parts is not None: + particles[popname].add(parts) + return particles + + +def amr_grid(hierarchy, time): + """returns a non-uniform contiguous primal grid + associated to the given hierarchy + """ + lvlPatchBoxes = {ilvl: [] for ilvl in range(hierarchy.finest_level() + 1)} + finalCells = {ilvl: None for ilvl in range(hierarchy.finest_level() + 1)} + lvl = hierarchy.levels(time) + + for ilvl in range(hierarchy.finest_level(time) + 1)[::-1]: + sorted_patches = sorted(lvl[ilvl].patches, key=lambda p: p.layout.box.lower[0]) + + for ip, patch in enumerate(sorted_patches): + box = patch.layout.box + lvlPatchBoxes[ilvl].append(box) + + # we create a list of all cells in the current patch + # remember that if the box upper cell is, say = 40, + # it means that the upper node is the lower node of cell 41 + # so to get all primal nodes of a patch we need to include + # one past the upper cell. + # this said we do not want to include that last primal nodes + # all the time because that would be a duplicate with the lower + # node of the next patch. We only want to add it for the LAST + # (because sorted) patch. We also do not want to do it on levels + # other than L0 because the last primal node of the last patch + # of L_i is the first primal node of a L_{i-1} node, so including it + # would also mean adding a duplicate. + last = 1 if ilvl == 0 and ip == len(sorted_patches) - 1 else 0 + cells = np.arange(box.lower[0], box.upper[0] + 1 + last) + + # finest level has no next finer so we take all cells + if ilvl == hierarchy.finest_level(time): + if finalCells[ilvl] is None: + finalCells[ilvl] = cells + else: + finalCells[ilvl] = np.concatenate((finalCells[ilvl], cells)) + + else: + # on other levels + # we take only grids not overlaped by next finer + coarsenedNextFinerBoxes = [ + boxm.coarsen(b, refinement_ratio) for b in lvlPatchBoxes[ilvl + 1] + ] + for coarseBox in coarsenedNextFinerBoxes: + ccells = np.arange(coarseBox.lower[0], coarseBox.upper[0] + 1) + inter, icells, iccells = np.intersect1d( + cells, ccells, return_indices=True + ) + cells = np.delete(cells, icells) + if len(cells): + if finalCells[ilvl] is None: + finalCells[ilvl] = cells + else: + finalCells[ilvl] = np.unique( + np.concatenate((finalCells[ilvl], cells)) + ) + + # now we have all cells for each level we + # just need to compute the primal coordinates + # and concatenate in a single array + for ilvl in range(hierarchy.finest_level() + 1): + if ilvl == 0: + x = finalCells[ilvl] * hierarchy.level(ilvl).patches[0].layout.dl[0] + else: + xx = finalCells[ilvl] * hierarchy.level(ilvl).patches[0].layout.dl[0] + x = np.concatenate((x, xx)) + + return np.sort(x) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py new file mode 100644 index 000000000..32c20722e --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py @@ -0,0 +1,402 @@ +from .hierarchy import PatchHierarchy +from .patchdata import FieldData +from .patchlevel import PatchLevel +from .patch import Patch +from ...core.phare_utilities import listify +from ...core.phare_utilities import refinement_ratio +from ...core import box as boxm + +import numpy as np + +field_qties = { + "EM_B_x": "Bx", + "EM_B_y": "By", + "EM_B_z": "Bz", + "EM_E_x": "Ex", + "EM_E_y": "Ey", + "EM_E_z": "Ez", + "flux_x": "Fx", + "flux_y": "Fy", + "flux_z": "Fz", + "bulkVelocity_x": "Vx", + "bulkVelocity_y": "Vy", + "bulkVelocity_z": "Vz", + "momentum_tensor_xx": "Mxx", + "momentum_tensor_xy": "Mxy", + "momentum_tensor_xz": "Mxz", + "momentum_tensor_yy": "Myy", + "momentum_tensor_yz": "Myz", + "momentum_tensor_zz": "Mzz", + "density": "rho", + "mass_density": "rho", + "tags": "tags", +} + + +def are_compatible_hierarchies(hierarchies): + return True + + +def merge_particles(hierarchy): + for time, patch_levels in hierarchy.time_hier.items(): + for ilvl, plvl in patch_levels.items(): + for ip, patch in enumerate(plvl.patches): + pdatas = patch.patch_datas + domain_pdata = [ + (pdname, pd) for pdname, pd in pdatas.items() if "domain" in pdname + ][0] + + pghost_pdatas = [ + (pdname, pd) + for pdname, pd in pdatas.items() + if "patchGhost" in pdname + ] + lghost_pdatas = [ + (pdname, pd) + for pdname, pd in pdatas.items() + if "levelGhost" in pdname + ] + + pghost_pdata = pghost_pdatas[0] if pghost_pdatas else None + lghost_pdata = lghost_pdatas[0] if lghost_pdatas else None + + if pghost_pdata is not None: + domain_pdata[1].dataset.add(pghost_pdata[1].dataset) + del pdatas[pghost_pdata[0]] + + if lghost_pdata is not None: + domain_pdata[1].dataset.add(lghost_pdata[1].dataset) + del pdatas[lghost_pdata[0]] + + popname = domain_pdata[0].split("_")[0] + pdatas[popname + "_particles"] = pdatas[domain_pdata[0]] + del pdatas[domain_pdata[0]] + + +def getPatch(hier, point): + """ + convenience function mainly for debug. returns a dict + {ilevel:patch} for patches in which the given point is + """ + patches = {} + counts = {ilvl: 0 for ilvl in hier.levels().keys()} + for ilvl, lvl in hier.levels().items(): + for p in lvl.patches: + px, py = point + dx, dy = p.layout.dl + ix = int(px / dx) + iy = int(py / dy) + if (ix, iy) in p.box: + patches[ilvl] = p + counts[ilvl] += 1 + for k, v in counts.items(): + if v > 1: + print("error : ", k, v) + raise RuntimeError("more than one patch found for point") + return patches + + +def compute_hier_from(compute, hierarchies, **kwargs): + """return a new hierarchy using the callback 'compute' on all patchdatas + of the given hierarchies + """ + if not are_compatible_hierarchies(hierarchies): + raise RuntimeError("hierarchies are not compatible") + + hierarchies = listify(hierarchies) + reference_hier = hierarchies[0] + domain_box = reference_hier.domain_box + patch_levels = {} + for ilvl in range(reference_hier.levelNbr()): + patch_levels[ilvl] = PatchLevel( + ilvl, new_patches_from(compute, hierarchies, ilvl, **kwargs) + ) + + assert len(reference_hier.time_hier) == 1 # only single time hierarchies now + t = list(reference_hier.time_hier.keys())[0] + return PatchHierarchy(patch_levels, domain_box, refinement_ratio, time=t) + + +def extract_patchdatas(hierarchies, ilvl, ipatch): + """ + returns a dict {patchdata_name:patchdata} from a list of hierarchies for patch ipatch at level ilvl + """ + patches = [h.patch_levels[ilvl].patches[ipatch] for h in hierarchies] + patch_datas = { + pdname: pd for p in patches for pdname, pd in list(p.patch_datas.items()) + } + return patch_datas + + +def new_patchdatas_from(compute, patchdatas, layout, id, **kwargs): + new_patch_datas = {} + datas = compute(patchdatas, id=id, **kwargs) + for data in datas: + pd = FieldData(layout, data["name"], data["data"], centering=data["centering"]) + new_patch_datas[data["name"]] = pd + return new_patch_datas + + +def new_patches_from(compute, hierarchies, ilvl, **kwargs): + reference_hier = hierarchies[0] + new_patches = [] + patch_nbr = len(reference_hier.patch_levels[ilvl].patches) + for ip in range(patch_nbr): + current_patch = reference_hier.patch_levels[ilvl].patches[ip] + layout = current_patch.layout + patch_datas = extract_patchdatas(hierarchies, ilvl, ip) + new_patch_datas = new_patchdatas_from( + compute, patch_datas, layout, id=current_patch.id, **kwargs + ) + new_patches.append(Patch(new_patch_datas, current_patch.id)) + return new_patches + + +def is_root_lvl(patch_level): + return patch_level.level_number == 0 + + +def quantidic(ilvl, wrangler): + pl = wrangler.getPatchLevel(ilvl) + + return { + "density": pl.getDensity, + "bulkVelocity_x": pl.getVix, + "bulkVelocity_y": pl.getViy, + "bulkVelocity_z": pl.getViz, + "EM_B_x": pl.getBx, + "EM_B_y": pl.getBy, + "EM_B_z": pl.getBz, + "EM_E_x": pl.getEx, + "EM_E_y": pl.getEy, + "EM_E_z": pl.getEz, + "flux_x": pl.getFx, + "flux_y": pl.getFy, + "flux_z": pl.getFz, + "particles": pl.getParticles, + } + + +def isFieldQty(qty): + return qty in ( + "density", + "bulkVelocity_x", + "bulkVelocity_y", + "bulkVelocity_z", + "EM_B_x", + "EM_B_y", + "EM_B_z", + "EM_E_x", + "EM_E_y", + "EM_E_z", + "flux_x", + "flux_y", + "flux_z", + "momentumTensor_xx", + "momentumTensor_xy", + "momentumTensor_xz", + "momentumTensor_yy", + "momentumTensor_yz", + "momentumTensor_zz", + ) + + +def overlap_mask_1d(x, dl, level, qty): + """ + return the mask for x where x is overlaped by the qty patch datas + on the given level, assuming that this level is finer than the one of x + + :param x: 1d array containing the [x] position + :param dl: list containing the grid steps where x is defined + :param level: a given level associated to a hierarchy + :param qty: ['Bx', 'By', 'Bz', 'Ex', 'Ey', 'Ez', 'Fx', 'Fy', 'Fz', 'Vx', 'Vy', 'Vz', 'rho'] + """ + + is_overlaped = np.ones(x.shape[0], dtype=bool) * False + + for patch in level.patches: + pdata = patch.patch_datas[qty] + ghosts_nbr = pdata.ghosts_nbr + + fine_x = pdata.x[ghosts_nbr[0] - 1 : -ghosts_nbr[0] + 1] + + fine_dl = pdata.dl + local_dl = dl + + if fine_dl[0] < local_dl[0]: + xmin, xmax = fine_x.min(), fine_x.max() + + overlaped_idx = np.where((x > xmin) & (x < xmax))[0] + + is_overlaped[overlaped_idx] = True + + else: + raise ValueError("level needs to have finer grid resolution than that of x") + + return is_overlaped + + +def overlap_mask_2d(x, y, dl, level, qty): + """ + return the mask for x & y where ix & y are overlaped by the qty patch datas + on the given level, assuming that this level is finer than the one of x & y + important note : this mask is flatten + + :param x: 1d array containing the [x] position + :param y: 1d array containing the [y] position + :param dl: list containing the grid steps where x and y are defined + :param level: a given level associated to a hierarchy + :param qty: ['Bx', 'By', 'Bz', 'Ex', 'Ey', 'Ez', 'Fx', 'Fy', 'Fz', 'Vx', 'Vy', 'Vz', 'rho'] + """ + + is_overlaped = np.ones([x.shape[0] * y.shape[0]], dtype=bool) * False + + for patch in level.patches: + pdata = patch.patch_datas[qty] + ghosts_nbr = pdata.ghosts_nbr + + fine_x = pdata.x[ghosts_nbr[0] - 1 : -ghosts_nbr[0] + 1] + fine_y = pdata.y[ghosts_nbr[1] - 1 : -ghosts_nbr[1] + 1] + + fine_dl = pdata.dl + local_dl = dl + + if (fine_dl[0] < local_dl[0]) and (fine_dl[1] < local_dl[1]): + xmin, xmax = fine_x.min(), fine_x.max() + ymin, ymax = fine_y.min(), fine_y.max() + + xv, yv = np.meshgrid(x, y, indexing="ij") + xf = xv.flatten() + yf = yv.flatten() + + overlaped_idx = np.where( + (xf > xmin) & (xf < xmax) & (yf > ymin) & (yf < ymax) + )[0] + + is_overlaped[overlaped_idx] = True + + else: + raise ValueError( + "level needs to have finer grid resolution than that of x or y" + ) + + return is_overlaped + + +def flat_finest_field(hierarchy, qty, time=None): + """ + returns 2 flattened arrays containing the data (with shape [Npoints]) + and the coordinates (with shape [Npoints, Ndim]) for the given + hierarchy of qty. + + :param hierarchy: the hierarchy where qty is defined + :param qty: the field (eg "Bx") that we want + """ + + dim = hierarchy.ndim + + if dim == 1: + return flat_finest_field_1d(hierarchy, qty, time) + elif dim == 2: + return flat_finest_field_2d(hierarchy, qty, time) + elif dim == 3: + raise RuntimeError("Not implemented") + # return flat_finest_field_3d(hierarchy, qty, time) + else: + raise ValueError("the dim of a hierarchy should be 1, 2 or 3") + + +def flat_finest_field_1d(hierarchy, qty, time=None): + lvl = hierarchy.levels(time) + + for ilvl in range(hierarchy.finest_level(time) + 1)[::-1]: + patches = lvl[ilvl].patches + + for ip, patch in enumerate(patches): + pdata = patch.patch_datas[qty] + + # all but 1 ghost nodes are removed in order to limit + # the overlapping, but to keep enough point to avoid + # any extrapolation for the interpolator + needed_points = pdata.ghosts_nbr - 1 + + # data = pdata.dataset[patch.box] # TODO : once PR 551 will be merged... + data = pdata.dataset[needed_points[0] : -needed_points[0]] + x = pdata.x[needed_points[0] : -needed_points[0]] + + if ilvl == hierarchy.finest_level(time): + if ip == 0: + final_data = data + final_x = x + else: + final_data = np.concatenate((final_data, data)) + final_x = np.concatenate((final_x, x)) + + else: + is_overlaped = overlap_mask_1d( + x, pdata.dl, hierarchy.level(ilvl + 1, time), qty + ) + + finest_data = data[~is_overlaped] + finest_x = x[~is_overlaped] + + final_data = np.concatenate((final_data, finest_data)) + final_x = np.concatenate((final_x, finest_x)) + + return final_data, final_x + + +def flat_finest_field_2d(hierarchy, qty, time=None): + lvl = hierarchy.levels(time) + + for ilvl in range(hierarchy.finest_level(time) + 1)[::-1]: + patches = lvl[ilvl].patches + + for ip, patch in enumerate(patches): + pdata = patch.patch_datas[qty] + + # all but 1 ghost nodes are removed in order to limit + # the overlapping, but to keep enough point to avoid + # any extrapolation for the interpolator + needed_points = pdata.ghosts_nbr - 1 + + # data = pdata.dataset[patch.box] # TODO : once PR 551 will be merged... + data = pdata.dataset[ + needed_points[0] : -needed_points[0], + needed_points[1] : -needed_points[1], + ] + x = pdata.x[needed_points[0] : -needed_points[0]] + y = pdata.y[needed_points[1] : -needed_points[1]] + + xv, yv = np.meshgrid(x, y, indexing="ij") + + data_f = data.flatten() + xv_f = xv.flatten() + yv_f = yv.flatten() + + if ilvl == hierarchy.finest_level(time): + if ip == 0: + final_data = data_f + tmp_x = xv_f + tmp_y = yv_f + else: + final_data = np.concatenate((final_data, data_f)) + tmp_x = np.concatenate((tmp_x, xv_f)) + tmp_y = np.concatenate((tmp_y, yv_f)) + + else: + is_overlaped = overlap_mask_2d( + x, y, pdata.dl, hierarchy.level(ilvl + 1, time), qty + ) + + finest_data = data_f[~is_overlaped] + finest_x = xv_f[~is_overlaped] + finest_y = yv_f[~is_overlaped] + + final_data = np.concatenate((final_data, finest_data)) + tmp_x = np.concatenate((tmp_x, finest_x)) + tmp_y = np.concatenate((tmp_y, finest_y)) + + final_xy = np.stack((tmp_x, tmp_y), axis=1) + + return final_data, final_xy diff --git a/pyphare/pyphare/pharesee/hierarchy/patch.py b/pyphare/pyphare/pharesee/hierarchy/patch.py new file mode 100644 index 000000000..75477ce68 --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/patch.py @@ -0,0 +1,66 @@ +from .patchdata import PatchData + + +class Patch: + """ + A patch represents a hyper-rectangular region of space + """ + + def __init__(self, patch_datas, patch_id="", layout=None, attrs=None): + """ + :param patch_datas: a list of PatchData objects + these are assumed to "belong" to the Patch so to + share the same origin, mesh size and box. + """ + if layout is not None: + self.layout = layout + self.box = layout.box + self.origin = layout.origin + self.dl = layout.dl + self.patch_datas = patch_datas + self.id = patch_id + + if len(patch_datas): + pdata0 = list(patch_datas.values())[0] # 0 represents all others + self.layout = pdata0.layout + self.box = pdata0.layout.box + self.origin = pdata0.layout.origin + self.dl = pdata0.layout.dl + self.patch_datas = patch_datas + self.id = patch_id + + self.attrs = attrs + + def __str__(self): + return f"Patch: box( {self.box}), id({self.id})" + + def __repr__(self): + return self.__str__() + + def copy(self): + """does not copy patchdatas.datasets (see class PatchData)""" + from copy import deepcopy + + return deepcopy(self) + + def __copy__(self): + return self.copy() + + def __call__(self, qty, **kwargs): + # take slice/slab of 1/2d array from 2/3d array + if "x" in kwargs and len(kwargs) == 1: + cut = kwargs["x"] + idim = 0 + elif "y" in kwargs and len(kwargs) == 1: + cut = kwargs["y"] + idim = 1 + else: + raise ValueError("need to specify either x or y cut coordinate") + pd = self.patch_datas[qty] + origin = pd.origin[idim] + idx = int((cut - origin) / pd.layout.dl[idim]) + nbrGhosts = pd.ghosts_nbr[idim] + if idim == 0: + return pd.dataset[idx + nbrGhosts, nbrGhosts:-nbrGhosts] + elif idim == 1: + return pd.dataset[nbrGhosts:-nbrGhosts, idx + nbrGhosts] diff --git a/pyphare/pyphare/pharesee/hierarchy/patchdata.py b/pyphare/pyphare/pharesee/hierarchy/patchdata.py new file mode 100644 index 000000000..22161d56e --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/patchdata.py @@ -0,0 +1,216 @@ +import numpy as np + +from ...core.phare_utilities import deep_copy +from ...core import box as boxm +from ...core.box import Box + + +class PatchData: + """ + base class for FieldData and ParticleData + this class just factors common geometrical properties + """ + + def __init__(self, layout, quantity): + """ + :param layout: a GridLayout representing the domain on which the data + is defined + :param quantity: ['field', 'particle'] + """ + self.quantity = quantity + self.box = layout.box + self.origin = layout.origin + self.layout = layout + + def __deepcopy__(self, memo): + no_copy_keys = ["dataset"] # do not copy these things + return deep_copy(self, memo, no_copy_keys) + + +class FieldData(PatchData): + """ + Concrete type of PatchData representing a physical quantity + defined on a grid. + """ + + @property + def x(self): + withGhosts = self.field_name != "tags" + if self._x is None: + self._x = self.layout.yeeCoordsFor( + self.field_name, + "x", + withGhosts=withGhosts, + centering=self.centerings[0], + ) + return self._x + + @property + def y(self): + withGhosts = self.field_name != "tags" + if self._y is None: + self._y = self.layout.yeeCoordsFor( + self.field_name, + "y", + withGhosts=withGhosts, + centering=self.centerings[1], + ) + return self._y + + @property + def z(self): + withGhosts = self.field_name != "tags" + if self._z is None: + self._z = self.layout.yeeCoordsFor( + self.field_name, + "z", + withGhosts=withGhosts, + centering=self.centerings[2], + ) + return self._z + + def primal_directions(self): + return self.size - self.ghost_box.shape + + def __str__(self): + return "FieldData: (box=({}, {}), key={})".format( + self.layout.box, self.layout.box.shape, self.field_name + ) + + def __repr__(self): + return self.__str__() + + def select(self, box): + """ + return view of internal data based on overlap of input box + returns a view +1 in size in primal directions + """ + assert isinstance(box, Box) and box.ndim == self.box.ndim + + gbox = self.ghost_box.copy() + gbox.upper += self.primal_directions() + + box = box.copy() + box.upper += self.primal_directions() + + overlap = box * gbox + if overlap is not None: + lower = self.layout.AMRToLocal(overlap.lower) + upper = self.layout.AMRToLocal(overlap.upper) + + if box.ndim == 1: + return self.dataset[lower[0] : upper[0] + 1] + if box.ndim == 2: + return self.dataset[lower[0] : upper[0] + 1, lower[1] : upper[1] + 1] + return np.array([]) + + def __getitem__(self, box): + return self.select(box) + + def __init__(self, layout, field_name, data, **kwargs): + """ + :param layout: A GridLayout representing the domain on which data is defined + :param field_name: the name of the field (e.g. "Bx") + :param data: the dataset from which data can be accessed + """ + super().__init__(layout, "field") + self._x = None + self._y = None + self._z = None + + self.layout = layout + self.field_name = field_name + self.name = field_name + self.dl = np.asarray(layout.dl) + self.ndim = layout.box.ndim + self.ghosts_nbr = np.zeros(self.ndim, dtype=int) + + if field_name in layout.centering["X"]: + directions = ["X", "Y", "Z"][: layout.box.ndim] # drop unused directions + self.centerings = [ + layout.qtyCentering(field_name, direction) for direction in directions + ] + elif "centering" in kwargs: + if isinstance(kwargs["centering"], list): + self.centerings = kwargs["centering"] + assert len(self.centerings) == self.ndim + else: + if self.ndim != 1: + raise ValueError( + "FieldData invalid dimenion for centering argument, expected list for dim > 1" + ) + self.centerings = [kwargs["centering"]] + else: + raise ValueError( + f"centering not specified and cannot be inferred from field name : {field_name}" + ) + + if self.field_name != "tags": + for i, centering in enumerate(self.centerings): + self.ghosts_nbr[i] = layout.nbrGhosts(layout.interp_order, centering) + + self.ghost_box = boxm.grow(layout.box, self.ghosts_nbr) + + self.size = np.copy(self.ghost_box.shape) + self.offset = np.zeros(self.ndim) + + for i, centering in enumerate(self.centerings): + if centering == "primal": + self.size[i] = self.ghost_box.shape[i] + 1 + else: + self.size[i] = self.ghost_box.shape[i] + self.offset[i] = 0.5 * self.dl[i] + + self.dataset = data + + def meshgrid(self, select=None): + def grid(): + if self.ndim == 1: + return [self.x] + if self.ndim == 2: + return np.meshgrid(self.x, self.y, indexing="ij") + return np.meshgrid(self.x, self.y, self.z, indexing="ij") + + mesh = grid() + if select is not None: + return tuple(g[select] for g in mesh) + return mesh + + +class ParticleData(PatchData): + """ + Concrete type of PatchData representing particles in a region + """ + + def __init__(self, layout, data, pop_name): + """ + :param layout: A GridLayout object representing the domain in which particles are + :param data: dataset containing particles + """ + super().__init__(layout, "particles") + self.dataset = data + self.pop_name = pop_name + self.name = pop_name + self.ndim = layout.box.ndim + + self.pop_name = pop_name + if layout.interp_order == 1: + self.ghosts_nbr = np.array([1] * layout.box.ndim) + elif layout.interp_order == 2 or layout.interp_order == 3: + self.ghosts_nbr = np.array([2] * layout.box.ndim) + else: + raise RuntimeError( + "invalid interpolation order {}".format(layout.interp_order) + ) + + self.ghost_box = boxm.grow(layout.box, self.ghosts_nbr) + assert (self.box.lower == self.ghost_box.lower + self.ghosts_nbr).all() + + def select(self, box): + return self.dataset[box] + + def __getitem__(self, box): + return self.select(box) + + def size(self): + return self.dataset.size() diff --git a/pyphare/pyphare/pharesee/hierarchy/patchlevel.py b/pyphare/pyphare/pharesee/hierarchy/patchlevel.py new file mode 100644 index 000000000..72fee0d5d --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/patchlevel.py @@ -0,0 +1,15 @@ +class PatchLevel: + """is a collection of patches""" + + def __init__(self, lvl_nbr, patches): + self.level_number = lvl_nbr + self.patches = patches + + def __iter__(self): + return self.patches.__iter__() + + def level_range(self): + name = list(self.patches[0].patch_datas.keys())[0] + return min([patch.patch_datas[name].x.min() for patch in self.patches]), max( + [patch.patch_datas[name].x.max() for patch in self.patches] + ) diff --git a/pyphare/pyphare/pharesee/plotting.py b/pyphare/pyphare/pharesee/plotting.py index 249d84624..9d73fa1bb 100644 --- a/pyphare/pyphare/pharesee/plotting.py +++ b/pyphare/pyphare/pharesee/plotting.py @@ -245,7 +245,7 @@ def finest_field_plot(run_path, qty, **kwargs): """ import os - from pyphare.pharesee.hierarchy import get_times_from_h5 + from pyphare.pharesee.hierarchy.fromh5 import get_times_from_h5 from pyphare.pharesee.run import Run from mpl_toolkits.axes_grid1 import make_axes_locatable import pyphare.core.gridlayout as gridlayout diff --git a/pyphare/pyphare/pharesee/run.py b/pyphare/pyphare/pharesee/run.py index 4fae4dedf..ac6d4ac69 100644 --- a/pyphare/pyphare/pharesee/run.py +++ b/pyphare/pyphare/pharesee/run.py @@ -1,11 +1,10 @@ import os import numpy as np -from .hierarchy import ( - compute_hier_from, - flat_finest_field, - hierarchy_from, -) +from .hierarchy.hierarchy_utils import compute_hier_from +from .hierarchy.hierarchy_utils import flat_finest_field +from .hierarchy import hierarchy_from + from pyphare.logger import getLogger logger = getLogger(__name__) diff --git a/pyphare/pyphare_tests/test_pharesee/__init__.py b/pyphare/pyphare_tests/test_pharesee/__init__.py index 5579bb6a6..621b7448f 100644 --- a/pyphare/pyphare_tests/test_pharesee/__init__.py +++ b/pyphare/pyphare_tests/test_pharesee/__init__.py @@ -7,10 +7,11 @@ from pyphare.pharesee.particles import Particles -from pyphare.pharesee.hierarchy import FieldData -from pyphare.pharesee.hierarchy import ParticleData +from pyphare.pharesee.hierarchy.patchdata import FieldData +from pyphare.pharesee.hierarchy.patchdata import ParticleData from pyphare.pharesee.hierarchy import PatchHierarchy -from pyphare.pharesee.hierarchy import Patch, PatchLevel +from pyphare.pharesee.hierarchy.patch import Patch +from pyphare.pharesee.hierarchy.patchlevel import PatchLevel """ number of ghosts is hard coded to 5 diff --git a/pyphare/pyphare_tests/test_pharesee/test_geometry.py b/pyphare/pyphare_tests/test_pharesee/test_geometry.py index a2153ce1b..c63a40b86 100644 --- a/pyphare/pyphare_tests/test_pharesee/test_geometry.py +++ b/pyphare/pyphare_tests/test_pharesee/test_geometry.py @@ -1,14 +1,7 @@ import unittest from ddt import ddt, data, unpack -import pyphare.core.box as boxm from pyphare.core.box import Box -from pyphare.core.phare_utilities import listify -from pyphare.pharesee.particles import Particles -from pyphare.pharesee.hierarchy import FieldData -from pyphare.pharesee.hierarchy import ParticleData -from pyphare.pharesee.hierarchy import PatchHierarchy -from pyphare.pharesee.hierarchy import Patch, PatchLevel from pyphare.pharesee.geometry import get_periodic_list, ghost_area_boxes from pyphare.pharesee.geometry import ( level_ghost_boxes, @@ -16,7 +9,7 @@ touch_domain_border, ) -from pyphare.core.gridlayout import GridLayout, yee_element_is_primal +from pyphare.core.gridlayout import yee_element_is_primal import numpy as np diff --git a/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py b/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py index 4dad87b1c..17af2eb2e 100644 --- a/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py +++ b/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py @@ -1,19 +1,11 @@ import unittest import numpy as np from ddt import ddt, data, unpack -import pyphare.core.box as boxm from pyphare.core.box import Box, Box2D -from pyphare.core.phare_utilities import listify -from pyphare.pharesee.particles import Particles -from pyphare.pharesee.hierarchy import FieldData -from pyphare.pharesee.hierarchy import ParticleData -from pyphare.pharesee.hierarchy import PatchHierarchy -from pyphare.pharesee.hierarchy import Patch, PatchLevel from pyphare.pharesee.geometry import ( level_ghost_boxes, hierarchy_overlaps, touch_domain_border, - toFieldBox, ghost_area_boxes, get_periodic_list, ) diff --git a/tests/amr/data/field/coarsening/test_coarsen_field.py b/tests/amr/data/field/coarsening/test_coarsen_field.py index 6068a3ba2..e9b3ed984 100644 --- a/tests/amr/data/field/coarsening/test_coarsen_field.py +++ b/tests/amr/data/field/coarsening/test_coarsen_field.py @@ -11,7 +11,7 @@ from pyphare.core.box import Box from pyphare.core.gridlayout import directions from pyphare.core.phare_utilities import refinement_ratio -from pyphare.pharesee.hierarchy import FieldData +from pyphare.pharesee.hierarchy.patchdata import FieldData def exec_fn(xyz, fn): diff --git a/tests/amr/data/field/refine/test_refine_field.py b/tests/amr/data/field/refine/test_refine_field.py index 88f0bbd77..188f80191 100644 --- a/tests/amr/data/field/refine/test_refine_field.py +++ b/tests/amr/data/field/refine/test_refine_field.py @@ -3,7 +3,7 @@ from pyphare.core import box as boxm from pyphare.core.gridlayout import GridLayout from pyphare.core.phare_utilities import refinement_ratio -from pyphare.pharesee.hierarchy import FieldData +from pyphare.pharesee.hierarchy.patchdata import FieldData # in this module, we assume the refinement ratio is 2 refinement_ratio = 2 diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py index 98701a09b..78e80bc38 100644 --- a/tests/simulator/test_advance.py +++ b/tests/simulator/test_advance.py @@ -17,7 +17,8 @@ ) from pyphare.pharein.simulation import Simulation from pyphare.pharesee.geometry import hierarchy_overlaps, level_ghost_boxes -from pyphare.pharesee.hierarchy import hierarchy_from, merge_particles +from pyphare.pharesee.hierarchy import hierarchy_from +from pyphare.pharesee.hierarchy.hierarchy_utils import merge_particles from pyphare.simulator.simulator import Simulator from tests.diagnostic import all_timestamps diff --git a/tests/simulator/test_diagnostics.py b/tests/simulator/test_diagnostics.py index 6b5ffe1d7..4814dc474 100644 --- a/tests/simulator/test_diagnostics.py +++ b/tests/simulator/test_diagnostics.py @@ -13,8 +13,9 @@ import pyphare.pharein as ph from ddt import data, ddt from pyphare.pharein.simulation import supported_dimensions -from pyphare.pharesee.hierarchy import h5_filename_from, h5_time_grp_key, hierarchy_from -from pyphare.simulator.simulator import Simulator, startMPI +from pyphare.pharesee.hierarchy.fromh5 import h5_filename_from, h5_time_grp_key +from pyphare.pharesee.hierarchy import hierarchy_from +from pyphare.simulator.simulator import Simulator from tests.diagnostic import dump_all_diags diff --git a/tests/simulator/test_initialization.py b/tests/simulator/test_initialization.py index 3604a3b89..7cfe0a960 100644 --- a/tests/simulator/test_initialization.py +++ b/tests/simulator/test_initialization.py @@ -16,7 +16,8 @@ ) from pyphare.pharein.simulation import Simulation from pyphare.pharesee.geometry import level_ghost_boxes -from pyphare.pharesee.hierarchy import hierarchy_from, merge_particles +from pyphare.pharesee.hierarchy.hierarchy_utils import merge_particles +from pyphare.pharesee.hierarchy import hierarchy_from from pyphare.pharesee.particles import aggregate as aggregate_particles from pyphare.simulator.simulator import Simulator @@ -718,7 +719,11 @@ def _test_patch_ghost_on_refined_level_case(self, dim, has_patch_ghost, **kwargs ] ) ) - self.assertTrue((1 in datahier.levels()) == has_patch_ghost) + nbrPatchGhostPatchDatasOnL1 = sum( + [len(p.patch_datas) for p in datahier.patch_levels[1].patches] + ) + + self.assertTrue((nbrPatchGhostPatchDatasOnL1 > 0) == has_patch_ghost) def _test_levelghostparticles_have_correct_split_from_coarser_particle( self, datahier From 0af8a15730439b67ebb3e28260ff3b821c8faea5 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Fri, 31 May 2024 13:57:18 +0200 Subject: [PATCH 03/18] split hierarchy_fromh5 (wip) --- pyphare/pyphare/pharesee/hierarchy/fromh5.py | 109 ++++++++++++++++-- .../pyphare/pharesee/hierarchy/hierarchy.py | 40 ++++++- 2 files changed, 135 insertions(+), 14 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/fromh5.py b/pyphare/pyphare/pharesee/hierarchy/fromh5.py index 710a82c2e..88f7a2313 100644 --- a/pyphare/pyphare/pharesee/hierarchy/fromh5.py +++ b/pyphare/pyphare/pharesee/hierarchy/fromh5.py @@ -10,6 +10,7 @@ from ...core.phare_utilities import refinement_ratio from ...core.gridlayout import GridLayout from .hierarchy_utils import field_qties +import h5py particle_files_patterns = ("domain", "patchGhost", "levelGhost") @@ -111,7 +112,6 @@ def h5_filename_from(diagInfo): def get_times_from_h5(filepath): - import h5py f = h5py.File(filepath, "r") times = np.array(sorted([float(s) for s in list(f["t"].keys())])) @@ -138,24 +138,109 @@ def load_one_time(time, hier): return time is not None and hier is not None -def fromh5_single_time(h5_filename, time, silent=True): - pass +def patch_levels_from_h5(filepath, time, selection_box=None): + """ + creates a dictionary of PatchLevels from a given time in a h5 file + {ilvl: PatchLevel} + """ + import os + h5f = h5py.File(filepath, "r") + root_cell_width = h5f.attrs["cell_width"] + interp_order = h5f.attrs["interpOrder"] + basename = os.path.basename(filepath) -def hierarchy_fromh5_(h5_filename, time, hier, silent=True): + patch_levels = {} - import h5py + for lvl_key, lvl in h5f["t"][time].items(): - data_file = h5py.File(h5_filename, "r") - basename = os.path.basename(h5_filename) - root_cell_width = np.asarray(data_file.attrs["cell_width"]) - interp = data_file.attrs["interpOrder"] - dimension = len(data_file.attrs["domain_box"]) - domain_box = Box([0] * dimension, data_file.attrs["domain_box"]) + ilvl = int(lvl_key[2:]) # pl1-->1 + lvl_cell_width = root_cell_width / refinement_ratio**ilvl + + patches = [] + + for h5_patch in lvl.values(): + lower = h5_patch.attrs["lower"] + upper = h5_patch.attrs["upper"] + origin = h5_patch.attrs["origin"] + + patch_box = Box(lower, upper) + + pos_upper = [ + orig + shape * dl + for orig, shape, dl in zip(origin, patch_box.shape, lvl_cell_width) + ] + pos_patch_box = Box(origin, pos_upper) + + intersect = None + if selection_box is not None: + intersect = selection_box * pos_patch_box + + if intersect is not None or selection_box is None: + patch_datas = {} + layout = make_layout(h5_patch, lvl_cell_width, interp_order) + if patch_has_datasets(h5_patch): + # we only add to patchdatas is there are datasets + # in the hdf5 patch group. + # but we do create a Patch (below) anyway since + # we want empty patches to be there as well to access their attributes. + add_to_patchdata(patch_datas, h5_patch, basename, layout) + + patches.append( + Patch( + patch_datas, + h5_patch.name.split("/")[-1], + layout=layout, + attrs={k: v for k, v in h5_patch.attrs.items()}, + ) + ) + + patch_levels[ilvl] = PatchLevel(ilvl, patches) + return patch_levels, h5f + + +def add_time_from_h5(hier, filepath, time, selection_box=None): + # add times to 'hier' + # we may have a different selection box for that time as for already existing times + # but we need to keep them, per time + if hier.has_time(time): + raise ValueError("time already exists in hierarchy") + + patch_levels, h5f = patch_levels_from_h5( + filepath, time, selection_box=selection_box + ) + + hier.add_time(time, patch_levels, h5f, selection_box=selection_box) + + return hier + + +def new_from_h5(filepath, time, selection_box=None): + # create a patchhierarchy from a given time and optional selection box + # loads all datasets from the filepath h5 file as patchdatas + + patch_levels, h5f = patch_levels_from_h5( + filepath, time, selection_box=selection_box + ) + dim = len(h5f.attrs["domain_box"]) + domain_box = Box([0] * dim, h5f.attrs["domain_box"]) + + # in hierarchy, we need to keep selection_box now + # because we want that operations involving several hierarchies will need to check + # that each time has the same patch layout. + hier = PatchHierarchy( + patch_levels, + domain_box, + refinement_ratio, + time, + h5f, + selection_box=selection_box, + ) + + return hier def hierarchy_fromh5(h5_filename, time, hier, silent=True): - import h5py data_file = h5py.File(h5_filename, "r") basename = os.path.basename(h5_filename) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index a1713dd3b..15702a780 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -3,6 +3,7 @@ from ...core.box import Box from ...core import box as boxm from ...core.phare_utilities import refinement_ratio +from ...core.phare_utilities import listify import numpy as np import matplotlib.pyplot as plt @@ -20,20 +21,40 @@ def __init__( data_files=None, **kwargs, ): + times = time + if not isinstance(times, (tuple, list)): + times = listify(times) + + if not isinstance(patch_levels, (tuple, list)): + patch_levels = listify(patch_levels) + + if "selection_box" in kwargs: + selection_box = kwargs["selection_box"] + if not isinstance(selection_box, (tuple, list)): + selection_box = listify(selection_box) + assert len(times) == len(selection_box) + + assert len(times) == len(patch_levels) + self.patch_levels = patch_levels self.ndim = len(domain_box.lower) self.time_hier = {} - self.time_hier.update({self.format_timestamp(time): patch_levels}) + self.time_hier.update( + {self.format_timestamp(t): pl for t, pl in zip(times, patch_levels)} + ) self.domain_box = domain_box self.refinement_ratio = refinement_ratio + self.selection_box = {self.format_timestamp(time): domain_box} self.data_files = {} self._sim = None if data_files is not None: - self.data_files.update(data_files) + self.data_files.update({data_files.filename: data_files}) + # TODO : need to update this after add_time() is called + # beter in another function if len(self.quantities()) > 1: for qty in self.quantities(): if qty in self.__dict__: @@ -170,6 +191,15 @@ def levelNbrs(self, time=None): time = self._default_time() return list(self.levels(time).keys()) + def add_time(self, time, patch_level, h5file, selection_box=None): + formated_time = self.format_timestamp(time) + + self.time_hier[self.format_timestamp(time)] = patch_level + if selection_box is not None: + self.selection_box[formated_time] = selection_box + + self.data_files[h5file.filename] = h5file + def is_homogeneous(self): """ return True if all patches of all levels at all times @@ -262,6 +292,12 @@ def __str__(self): s = s + "\n" return s + def has_time(self, time): + return self.format_timestamp(time) in self.time_hier + + def has_file(self, filename): + return filename in self.data_files + def times(self): return np.sort(np.asarray(list(self.time_hier.keys()))) From c5c75316d790fa7da18b06babbedc1a357fa594b Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Fri, 31 May 2024 18:47:02 +0200 Subject: [PATCH 04/18] fix data_files and particle count phloping --- .../pyphare/pharesee/hierarchy/hierarchy.py | 8 ++++-- .../test_pharesee/test_geometry_2d.py | 8 +++--- tools/python3/phloping.py | 28 ++++++++----------- 3 files changed, 22 insertions(+), 22 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index 15702a780..be5f33dec 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -83,14 +83,18 @@ def sim(self): if self._sim: return self._sim - if "py_attrs" not in self.data_files: + # data_files has a key/value per h5 filename. + # but the "serialized_simulation" in "py_attrs" should be the same for all files + # used by the hierarchy. So we just take the first one. + first_file = list(self.data_files.values())[0] + if "py_attrs" not in first_file.keys(): raise ValueError("Simulation is not available for deserialization") from ...pharein.simulation import deserialize try: self._sim = deserialize( - self.data_files["py_attrs"].attrs["serialized_simulation"] + first_file["py_attrs"].attrs["serialized_simulation"] ) except Exception as e: raise RuntimeError(f"Failed to deserialize simulation from data file : {e}") diff --git a/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py b/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py index 17af2eb2e..48eb1bcc1 100644 --- a/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py +++ b/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py @@ -115,7 +115,7 @@ def test_overlaps(self, refinement_boxes, expected): ) level_overlaps = hierarchy_overlaps(hierarchy) - for ilvl, lvl in enumerate(hierarchy.patch_levels): + for ilvl, lvl in enumerate(hierarchy.levels().items()): if ilvl not in expected: continue self.assertEqual(len(expected[ilvl]), len(level_overlaps[ilvl])) @@ -267,7 +267,7 @@ def test_particle_ghost_area_boxes(self): self.assertEqual(len(expected), len(gaboxes)) - for ilvl, lvl in enumerate(hierarchy.patch_levels): + for ilvl, lvl in enumerate(hierarchy.levels().items()): self.assertEqual(len(gaboxes[ilvl][particles]), len(expected[ilvl])) for act_pdata, exp_pdata in zip(gaboxes[ilvl][particles], expected[ilvl]): self.assertEqual(len(exp_pdata["boxes"]), len(act_pdata["boxes"])) @@ -289,7 +289,7 @@ def test_particle_level_ghost_boxes_do_not_overlap_patch_interiors(self): assert len(gaboxes_list) > 0 for pdatainfo in gaboxes_list: for box in pdatainfo["boxes"]: - for patch in hierarchy.patch_levels[ilvl].patches: + for patch in hierarchy.level(ilvl).patches: self.assertIsNone(patch.box * box) @@ -378,7 +378,7 @@ def test_level_ghost_boxes(self, refinement_boxes, expected): ) lvl_gaboxes = level_ghost_boxes(hierarchy, "particles") - for ilvl in range(1, len(hierarchy.patch_levels)): + for ilvl in range(1, len(hierarchy.levels())): self.assertEqual(len(lvl_gaboxes[ilvl].keys()), 1) key = list(lvl_gaboxes[ilvl].keys())[0] diff --git a/tools/python3/phloping.py b/tools/python3/phloping.py index fd980dfd8..99eaabee2 100644 --- a/tools/python3/phloping.py +++ b/tools/python3/phloping.py @@ -38,25 +38,21 @@ def _particles(self): """ Extract particle count per timestep per level from phlop logging """ + from pyphare.pharesee.hierarchy.fromh5 import get_times_from_h5 + bad_diag_error = ( "Simulation was not configured with Particle Count info diagnostic" ) - pcount_hier = hierarchy_from(h5_filename=self.run.path + "/particle_count.h5") - particles_per_level_per_time_step = { # per coarse timestep only - li: np.zeros(len(self.advances)) - for li in range(len(pcount_hier.data_files["t"][pcount_hier.times()[0]])) - } - for ti, t in enumerate(pcount_hier.times()[1:]): - for plk, pl in pcount_hier.data_files["t"][t].items(): - pc = 0 - for pid, p in pl.items(): - if "particle_count" not in p.attrs: - raise ValueError(bad_diag_error) - if pid.split("#")[0][1:] == self.rank: - pc += p.attrs["particle_count"] - if pc == 0: - pc += 1 # avoid div by 0 for rank missing patch on level - particles_per_level_per_time_step[int(plk[2:])][ti] += pc + filepath = self.run.path + "/particle_count.h5" + all_times = get_times_from_h5(filepath) + + particles_per_level_per_time_step = {} + pcount_hier = None + for it, time in enumerate(all_times): + pcount_hier = self.run.GetParticleCount(time) + for ilvl, lvl in pcount_hier.levels().items(): + pc = sum([p.attrs["particle_count"] for p in lvl.patches]) + particles_per_level_per_time_step[ilvl][it] = pc return pcount_hier, particles_per_level_per_time_step def advance_times_for_L(self, ilvl): From ea5e89841ff09820c5340b92ca0845fe4958b85b Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Fri, 31 May 2024 22:26:42 +0200 Subject: [PATCH 05/18] add new data to existing hierarchy --- pyphare/pyphare/pharesee/hierarchy/fromh5.py | 22 +++++++++++++++++++ .../pyphare/pharesee/hierarchy/hierarchy.py | 4 +++- tools/python3/phloping.py | 2 +- 3 files changed, 26 insertions(+), 2 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/fromh5.py b/pyphare/pyphare/pharesee/hierarchy/fromh5.py index 88f7a2313..b215eda51 100644 --- a/pyphare/pyphare/pharesee/hierarchy/fromh5.py +++ b/pyphare/pyphare/pharesee/hierarchy/fromh5.py @@ -215,6 +215,28 @@ def add_time_from_h5(hier, filepath, time, selection_box=None): return hier +def add_data_from_h5(hier, filepath, time): + """ + adds new PatchDatas to an existing PatchHierarchy for an existing time + + Data will be added from the given filepath. + Data will be extracted from the selection box of the hierarchy at that time. + """ + if not hier.has_time(time): + raise ValueError("time does not exist in hierarchy") + + # force using the hierarchy selection box at that time if existing + patch_levels, h5f = patch_levels_from_h5( + filepath, time, selection_box=hier.selection_box[time] + ) + + for ilvl, lvl in hier.levels(time).items(): + for ip, patch in enumerate(lvl.patches): + patch.patch_datas.update(patch_levels[ilvl].patches[ip].patch_datas) + + return hier + + def new_from_h5(filepath, time, selection_box=None): # create a patchhierarchy from a given time and optional selection box # loads all datasets from the filepath h5 file as patchdatas diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index be5f33dec..1821390f4 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -45,7 +45,9 @@ def __init__( self.domain_box = domain_box self.refinement_ratio = refinement_ratio - self.selection_box = {self.format_timestamp(time): domain_box} + self.selection_box = { + self.format_timestamp(t): box for t, box in zip(times, selection_box) + } self.data_files = {} self._sim = None diff --git a/tools/python3/phloping.py b/tools/python3/phloping.py index 99eaabee2..ac74beb26 100644 --- a/tools/python3/phloping.py +++ b/tools/python3/phloping.py @@ -50,7 +50,7 @@ def _particles(self): pcount_hier = None for it, time in enumerate(all_times): pcount_hier = self.run.GetParticleCount(time) - for ilvl, lvl in pcount_hier.levels().items(): + for ilvl, lvl in pcount_hier.levels(time).items(): pc = sum([p.attrs["particle_count"] for p in lvl.patches]) particles_per_level_per_time_step[ilvl][it] = pc return pcount_hier, particles_per_level_per_time_step From a736996828400380cbeb357926a384eb3554ce80 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Fri, 31 May 2024 22:45:04 +0200 Subject: [PATCH 06/18] fix misc issues --- pyphare/pyphare/pharesee/hierarchy/hierarchy.py | 6 +++--- tests/functional/alfven_wave/alfven_wave1d.py | 2 +- tests/functional/conservation/conserv.py | 12 +++++------- tests/simulator/test_advance.py | 3 ++- tests/simulator/test_initialization.py | 2 +- tests/simulator/test_restarts.py | 4 ++-- 6 files changed, 14 insertions(+), 15 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index 1821390f4..b53ff81da 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -32,6 +32,9 @@ def __init__( selection_box = kwargs["selection_box"] if not isinstance(selection_box, (tuple, list)): selection_box = listify(selection_box) + self.selection_box = { + self.format_timestamp(t): box for t, box in zip(times, selection_box) + } assert len(times) == len(selection_box) assert len(times) == len(patch_levels) @@ -45,9 +48,6 @@ def __init__( self.domain_box = domain_box self.refinement_ratio = refinement_ratio - self.selection_box = { - self.format_timestamp(t): box for t, box in zip(times, selection_box) - } self.data_files = {} self._sim = None diff --git a/tests/functional/alfven_wave/alfven_wave1d.py b/tests/functional/alfven_wave/alfven_wave1d.py index aa9ee6846..c9d72909a 100644 --- a/tests/functional/alfven_wave/alfven_wave1d.py +++ b/tests/functional/alfven_wave/alfven_wave1d.py @@ -161,7 +161,7 @@ def main(): by, xby = flat_finest_field(B, "By") ax.plot(xby, by, label="t = 500", alpha=0.6) - sorted_patches = sorted(B.patch_levels[1].patches, key=lambda p: p.box.lower[0]) + sorted_patches = sorted(B.level(1).patches, key=lambda p: p.box.lower[0]) x0 = sorted_patches[0].patch_datas["By"].x[0] x1 = sorted_patches[-1].patch_datas["By"].x[-1] diff --git a/tests/functional/conservation/conserv.py b/tests/functional/conservation/conserv.py index e0c5773c8..6e30b3dd1 100644 --- a/tests/functional/conservation/conserv.py +++ b/tests/functional/conservation/conserv.py @@ -131,10 +131,10 @@ def total_particles(parts, fun, lvlNbr=0, **kwargs): quantity on a given level, quantity being estimated by the callback "fun" (could be kinetic energy, momentum, et.) """ - for ilvl, lvl in parts.patch_levels.items(): + tot = 0.0 + for ilvl, lvl in parts.levels().items(): if lvlNbr == ilvl: - tot = 0.0 - for ip, patch in enumerate(lvl.patches): + for patch in lvl.patches: keys = list(patch.patch_datas.keys()) pdata = patch.patch_datas[keys[0]] particles = pdata.dataset @@ -151,7 +151,7 @@ def mag_energy(B, lvlNbr=0): """ return the total magnetic energy on a given level """ - for ilvl, lvl in B.patch_levels.items(): + for ilvl, lvl in B.levels().items(): if lvlNbr == ilvl: tot = 0.0 for ip, patch in enumerate(lvl.patches): @@ -172,9 +172,7 @@ def mag_energy(B, lvlNbr=0): bz = 0.5 * (bztmp[1:] + bztmp[:-1]) # sum 0.5B^2 * dx over all nodes - per_patch = np.sum( - (bx**2 + by**2 + bz**2) * 0.5 * pdata.layout.dl[0] - ) + per_patch = np.sum((bx**2 + by**2 + bz**2) * 0.5 * pdata.layout.dl[0]) tot += per_patch return tot diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py index 78e80bc38..4d07aa147 100644 --- a/tests/simulator/test_advance.py +++ b/tests/simulator/test_advance.py @@ -357,7 +357,8 @@ def _test_overlapped_particledatas_have_identical_particles( overlaps = hierarchy_overlaps(datahier, coarsest_time) - for ilvl, lvl in datahier.patch_levels.items(): + # should it be levels(coarsest_time)? + for ilvl, lvl in datahier.levels().items(): print("testing level {}".format(ilvl)) for overlap in overlaps[ilvl]: pd1, pd2 = overlap["pdatas"] diff --git a/tests/simulator/test_initialization.py b/tests/simulator/test_initialization.py index 7cfe0a960..9f4ed8ec4 100644 --- a/tests/simulator/test_initialization.py +++ b/tests/simulator/test_initialization.py @@ -720,7 +720,7 @@ def _test_patch_ghost_on_refined_level_case(self, dim, has_patch_ghost, **kwargs ) ) nbrPatchGhostPatchDatasOnL1 = sum( - [len(p.patch_datas) for p in datahier.patch_levels[1].patches] + [len(p.patch_datas) for p in datahier.level(1).patches] ) self.assertTrue((nbrPatchGhostPatchDatasOnL1 > 0) == has_patch_ghost) diff --git a/tests/simulator/test_restarts.py b/tests/simulator/test_restarts.py index f847ac56c..e46b5af25 100644 --- a/tests/simulator/test_restarts.py +++ b/tests/simulator/test_restarts.py @@ -155,8 +155,8 @@ def check_diags(self, diag_dir0, diag_dir1, pops, timestamps, expected_num_level def count_levels_and_patches(qty): n_levels = len(qty.levels()) n_patches = 0 - for ilvl, lvl in qty.patch_levels.items(): - n_patches += len(qty.patch_levels[ilvl].patches) + for ilvl in qty.levels().keys(): + n_patches += len(qty.level(ilvl).patches) return n_levels, n_patches self.assertGreater(len(timestamps), 0) From 281e4135cb84ffecfb9cca0301cfd1c1e54a0b85 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Sat, 1 Jun 2024 13:02:17 +0200 Subject: [PATCH 07/18] multiple times compute_hier_from --- .../pharesee/hierarchy/hierarchy_utils.py | 45 ++++++++++++------- 1 file changed, 29 insertions(+), 16 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py index 32c20722e..e63093915 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py @@ -34,6 +34,16 @@ def are_compatible_hierarchies(hierarchies): + ref = hierarchies[0] + same_box = True + same_selection = True + same_files = True + same_times = True + for hier in hierarchies[1:]: + same_box = same_box and (hier.domain_box == ref.domain_box) + same_selection = same_selection and (hier.selection_box == ref.selection_box) + same_files = same_files and (hier.data_files.keys() == ref.data_files.keys()) + same_times = same_times and (hier.time_hier.keys() == ref.time_hier.keys()) return True @@ -100,28 +110,31 @@ def compute_hier_from(compute, hierarchies, **kwargs): """return a new hierarchy using the callback 'compute' on all patchdatas of the given hierarchies """ - if not are_compatible_hierarchies(hierarchies): - raise RuntimeError("hierarchies are not compatible") hierarchies = listify(hierarchies) + if not are_compatible_hierarchies(hierarchies): + raise RuntimeError("hierarchies are not compatible") reference_hier = hierarchies[0] domain_box = reference_hier.domain_box - patch_levels = {} - for ilvl in range(reference_hier.levelNbr()): - patch_levels[ilvl] = PatchLevel( - ilvl, new_patches_from(compute, hierarchies, ilvl, **kwargs) - ) + patch_levels_per_time = [] + for t in reference_hier.times(): + patch_levels = {} + for ilvl in range(reference_hier.levelNbr()): + patch_levels[ilvl] = PatchLevel( + ilvl, new_patches_from(compute, hierarchies, ilvl, t, **kwargs) + ) + patch_levels_per_time.append(patch_levels) - assert len(reference_hier.time_hier) == 1 # only single time hierarchies now - t = list(reference_hier.time_hier.keys())[0] - return PatchHierarchy(patch_levels, domain_box, refinement_ratio, time=t) + return PatchHierarchy( + patch_levels_per_time, domain_box, refinement_ratio, time=reference_hier.times() + ) -def extract_patchdatas(hierarchies, ilvl, ipatch): +def extract_patchdatas(hierarchies, ilvl, t, ipatch): """ returns a dict {patchdata_name:patchdata} from a list of hierarchies for patch ipatch at level ilvl """ - patches = [h.patch_levels[ilvl].patches[ipatch] for h in hierarchies] + patches = [h.level(ilvl, time=t).patches[ipatch] for h in hierarchies] patch_datas = { pdname: pd for p in patches for pdname, pd in list(p.patch_datas.items()) } @@ -137,14 +150,14 @@ def new_patchdatas_from(compute, patchdatas, layout, id, **kwargs): return new_patch_datas -def new_patches_from(compute, hierarchies, ilvl, **kwargs): +def new_patches_from(compute, hierarchies, ilvl, t, **kwargs): reference_hier = hierarchies[0] new_patches = [] - patch_nbr = len(reference_hier.patch_levels[ilvl].patches) + patch_nbr = len(reference_hier.level(ilvl, time=t).patches) for ip in range(patch_nbr): - current_patch = reference_hier.patch_levels[ilvl].patches[ip] + current_patch = reference_hier.level(ilvl, time=t).patches[ip] layout = current_patch.layout - patch_datas = extract_patchdatas(hierarchies, ilvl, ip) + patch_datas = extract_patchdatas(hierarchies, ilvl, t, ip) new_patch_datas = new_patchdatas_from( compute, patch_datas, layout, id=current_patch.id, **kwargs ) From a60e10a74ee34a03c99d3101ed59fbba700d0653 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Sat, 1 Jun 2024 14:34:49 +0200 Subject: [PATCH 08/18] up --- pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py index e63093915..d61989561 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py @@ -4,7 +4,6 @@ from .patch import Patch from ...core.phare_utilities import listify from ...core.phare_utilities import refinement_ratio -from ...core import box as boxm import numpy as np From c26dde2c99befd2cd657ef169085ce2fa0336d77 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Sat, 1 Jun 2024 14:46:34 +0200 Subject: [PATCH 09/18] update hierarchy attributes --- pyphare/pyphare/pharesee/hierarchy/hierarchy.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index b53ff81da..c841b2ffb 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -55,8 +55,12 @@ def __init__( if data_files is not None: self.data_files.update({data_files.filename: data_files}) - # TODO : need to update this after add_time() is called - # beter in another function + self.update() + + def __getitem__(self, qty): + return self.__dict__[qty] + + def update(self): if len(self.quantities()) > 1: for qty in self.quantities(): if qty in self.__dict__: @@ -77,9 +81,6 @@ def __init__( else: self.qty.time_hier[time] = new_lvls # pylint: disable=E1101 - def __getitem__(self, qty): - return self.__dict__[qty] - @property def sim(self): if self._sim: @@ -205,6 +206,7 @@ def add_time(self, time, patch_level, h5file, selection_box=None): self.selection_box[formated_time] = selection_box self.data_files[h5file.filename] = h5file + self.update() def is_homogeneous(self): """ From 875361473d04e8c360ecb026994e33809bf3220c Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Sat, 1 Jun 2024 16:38:55 +0200 Subject: [PATCH 10/18] wip HFH5 --- pyphare/pyphare/pharesee/hierarchy/fromh5.py | 21 ++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/pyphare/pyphare/pharesee/hierarchy/fromh5.py b/pyphare/pyphare/pharesee/hierarchy/fromh5.py index b215eda51..d9f7ea886 100644 --- a/pyphare/pyphare/pharesee/hierarchy/fromh5.py +++ b/pyphare/pyphare/pharesee/hierarchy/fromh5.py @@ -262,6 +262,27 @@ def new_from_h5(filepath, time, selection_box=None): return hier +def hierarchy_fromh5_(h5_filename, time=None, hier=None, silent=True): + """ + creates a PatchHierarchy from a given time in a h5 file + if hier is None, a new hierarchy is created + if hier is not None, data is added to the hierarchy + """ + if create_from_one_time(time, hier): + return new_from_h5(h5_filename, time) + + if create_from_all_times(time, hier): + times = get_times_from_h5(h5_filename) + h = new_from_h5(h5_filename, times[0]) + for time in times[1:]: + add_time_from_h5(h, h5_filename, time) + + if load_one_time(time, hier): + return add_time_from_h5(hier, h5_filename, time) + + return add_data_from_h5(hier, h5_filename, time) + + def hierarchy_fromh5(h5_filename, time, hier, silent=True): data_file = h5py.File(h5_filename, "r") From feb31e314ec0cb0fd07fc9ff500a9711c655d031 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Sat, 1 Jun 2024 17:33:18 +0200 Subject: [PATCH 11/18] update hierarchy quantity members --- pyphare/pyphare/pharesee/hierarchy/hierarchy.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index c841b2ffb..53362d178 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -63,9 +63,6 @@ def __getitem__(self, qty): def update(self): if len(self.quantities()) > 1: for qty in self.quantities(): - if qty in self.__dict__: - continue - first = True for time, levels in self.time_hier.items(): new_lvls = {} for ilvl, level in levels.items(): @@ -73,13 +70,12 @@ def update(self): for patch in level.patches: patches += [Patch({qty: patch.patch_datas[qty]}, patch.id)] new_lvls[ilvl] = PatchLevel(ilvl, patches) - if first: + if qty not in self.__dict__: self.__dict__[qty] = PatchHierarchy( new_lvls, self.domain_box, time=time ) - first = False else: - self.qty.time_hier[time] = new_lvls # pylint: disable=E1101 + self.__dict__[qty].time_hier[time] = new_lvls @property def sim(self): From 640c5e3768774169d9fb6fdfc67702d5f8581ee8 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Sun, 2 Jun 2024 18:00:35 +0200 Subject: [PATCH 12/18] up --- pyphare/pyphare/pharesee/hierarchy/fromh5.py | 2 +- .../pyphare/pharesee/hierarchy/hierarchy.py | 29 ++++++++++++------- 2 files changed, 20 insertions(+), 11 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/fromh5.py b/pyphare/pyphare/pharesee/hierarchy/fromh5.py index d9f7ea886..73d3027f7 100644 --- a/pyphare/pyphare/pharesee/hierarchy/fromh5.py +++ b/pyphare/pyphare/pharesee/hierarchy/fromh5.py @@ -343,7 +343,7 @@ def hierarchy_fromh5(h5_filename, time, hier, silent=True): patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl]) diag_hier = PatchHierarchy( - patch_levels, domain_box, refinement_ratio, t, data_file + patch_levels, domain_box, refinement_ratio, t, data_files=data_file ) return diag_hier diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index 53362d178..f2044761b 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -28,14 +28,15 @@ def __init__( if not isinstance(patch_levels, (tuple, list)): patch_levels = listify(patch_levels) - if "selection_box" in kwargs: - selection_box = kwargs["selection_box"] - if not isinstance(selection_box, (tuple, list)): - selection_box = listify(selection_box) + self.selection_box = kwargs.get("selection_box", None) + if self.selection_box is not None: + if not isinstance(self.selection_box, (tuple, list)): + self.selection_box = listify(self.selection_box) self.selection_box = { - self.format_timestamp(t): box for t, box in zip(times, selection_box) + self.format_timestamp(t): box + for t, box in zip(times, self.selection_box) } - assert len(times) == len(selection_box) + assert len(times) == len(self.selection_box) assert len(times) == len(patch_levels) @@ -49,11 +50,15 @@ def __init__( self.domain_box = domain_box self.refinement_ratio = refinement_ratio - self.data_files = {} self._sim = None - if data_files is not None: - self.data_files.update({data_files.filename: data_files}) + if data_files is not None and isinstance(data_files, dict): + self.data_files = data_files + elif data_files is not None: + if hasattr(self, "data_files"): + self.data_files.update({data_files.filename: data_files}) + else: + self.data_files = {data_files.filename: data_files} self.update() @@ -72,7 +77,11 @@ def update(self): new_lvls[ilvl] = PatchLevel(ilvl, patches) if qty not in self.__dict__: self.__dict__[qty] = PatchHierarchy( - new_lvls, self.domain_box, time=time + new_lvls, + self.domain_box, + selection_box=self.domain_box, + time=time, + data_files=self.data_files, ) else: self.__dict__[qty].time_hier[time] = new_lvls From cc92466de17e6181e327cf81b40a1a5b796f7b06 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Tue, 4 Jun 2024 11:43:22 +0200 Subject: [PATCH 13/18] up --- .../pyphare/pharesee/hierarchy/__init__.py | 4 +- pyphare/pyphare/pharesee/hierarchy/fromh5.py | 47 ++++---- .../pyphare/pharesee/hierarchy/hierarchy.py | 19 ++- pyphare/pyphare/pharesee/run.py | 108 ++++++++++-------- 4 files changed, 104 insertions(+), 74 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/__init__.py b/pyphare/pyphare/pharesee/hierarchy/__init__.py index edeed0416..c892caec3 100644 --- a/pyphare/pyphare/pharesee/hierarchy/__init__.py +++ b/pyphare/pyphare/pharesee/hierarchy/__init__.py @@ -4,7 +4,7 @@ def hierarchy_from( - simulator=None, qty=None, pop="", h5_filename=None, time=None, hier=None + simulator=None, qty=None, pop="", h5_filename=None, times=None, hier=None, **kwargs ): from .fromh5 import hierarchy_fromh5 from .fromsim import hierarchy_from_sim @@ -24,7 +24,7 @@ def hierarchy_from( raise ValueError("cannot pass both a simulator and a h5 file") if h5_filename is not None: - return hierarchy_fromh5(h5_filename, time, hier) + return hierarchy_fromh5(h5_filename, times, hier, **kwargs) if simulator is not None and qty is not None: return hierarchy_from_sim(simulator, qty, pop=pop) diff --git a/pyphare/pyphare/pharesee/hierarchy/fromh5.py b/pyphare/pyphare/pharesee/hierarchy/fromh5.py index 73d3027f7..a8be6fd1f 100644 --- a/pyphare/pyphare/pharesee/hierarchy/fromh5.py +++ b/pyphare/pyphare/pharesee/hierarchy/fromh5.py @@ -126,8 +126,8 @@ def create_from_all_times(time, hier): return time is None and hier is None -def create_from_one_time(time, hier): - return time is not None and hier is None +def create_from_times(times, hier): + return times is not None and hier is None def load_all_times(time, hier): @@ -138,17 +138,16 @@ def load_one_time(time, hier): return time is not None and hier is not None -def patch_levels_from_h5(filepath, time, selection_box=None): +def patch_levels_from_h5(h5f, time, selection_box=None): """ creates a dictionary of PatchLevels from a given time in a h5 file {ilvl: PatchLevel} """ import os - h5f = h5py.File(filepath, "r") root_cell_width = h5f.attrs["cell_width"] interp_order = h5f.attrs["interpOrder"] - basename = os.path.basename(filepath) + basename = os.path.basename(h5f.filename) patch_levels = {} @@ -196,13 +195,15 @@ def patch_levels_from_h5(filepath, time, selection_box=None): ) patch_levels[ilvl] = PatchLevel(ilvl, patches) - return patch_levels, h5f + return patch_levels -def add_time_from_h5(hier, filepath, time, selection_box=None): +def add_time_from_h5(hier, filepath, time, **kwargs): # add times to 'hier' # we may have a different selection box for that time as for already existing times # but we need to keep them, per time + + selection_box = kwargs.get("selection_box", None) if hier.has_time(time): raise ValueError("time already exists in hierarchy") @@ -237,13 +238,19 @@ def add_data_from_h5(hier, filepath, time): return hier -def new_from_h5(filepath, time, selection_box=None): +def new_from_h5(filepath, times, **kwargs): # create a patchhierarchy from a given time and optional selection box # loads all datasets from the filepath h5 file as patchdatas - patch_levels, h5f = patch_levels_from_h5( - filepath, time, selection_box=selection_box - ) + selection_box = kwargs.get("selection_box", None) + + patch_levels_per_time = [] + + h5f = h5py.File(filepath, "r") + for time in times: + patch_levels = patch_levels_from_h5(h5f, time, selection_box=selection_box) + patch_levels_per_time.append(patch_levels) + dim = len(h5f.attrs["domain_box"]) domain_box = Box([0] * dim, h5f.attrs["domain_box"]) @@ -251,10 +258,10 @@ def new_from_h5(filepath, time, selection_box=None): # because we want that operations involving several hierarchies will need to check # that each time has the same patch layout. hier = PatchHierarchy( - patch_levels, + patch_levels_per_time, domain_box, refinement_ratio, - time, + times, h5f, selection_box=selection_box, ) @@ -262,28 +269,28 @@ def new_from_h5(filepath, time, selection_box=None): return hier -def hierarchy_fromh5_(h5_filename, time=None, hier=None, silent=True): +def hierarchy_fromh5(h5_filename, time=None, hier=None, silent=True, **kwargs): """ creates a PatchHierarchy from a given time in a h5 file if hier is None, a new hierarchy is created if hier is not None, data is added to the hierarchy """ - if create_from_one_time(time, hier): - return new_from_h5(h5_filename, time) + if create_from_times(time, hier): + return new_from_h5(h5_filename, time, **kwargs) if create_from_all_times(time, hier): times = get_times_from_h5(h5_filename) - h = new_from_h5(h5_filename, times[0]) + h = new_from_h5(h5_filename, times[0], **kwargs) for time in times[1:]: - add_time_from_h5(h, h5_filename, time) + add_time_from_h5(h, h5_filename, time, **kwargs) if load_one_time(time, hier): - return add_time_from_h5(hier, h5_filename, time) + return add_time_from_h5(hier, h5_filename, time, **kwargs) return add_data_from_h5(hier, h5_filename, time) -def hierarchy_fromh5(h5_filename, time, hier, silent=True): +def hierarchy_fromh5_(h5_filename, time, hier, silent=True): data_file = h5py.File(h5_filename, "r") basename = os.path.basename(h5_filename) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index f2044761b..ba3b6fd48 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -17,11 +17,10 @@ def __init__( patch_levels, domain_box, refinement_ratio=2, - time=0.0, + times=[0.0], data_files=None, **kwargs, ): - times = time if not isinstance(times, (tuple, list)): times = listify(times) @@ -86,6 +85,22 @@ def update(self): else: self.__dict__[qty].time_hier[time] = new_lvls + def nbytes(self): + n = 0 + for t in self.times(): + for lvl in self.levels(t).values(): + for p in lvl.patches: + for pd in p.patch_datas.values(): + n += pd.dataset.nbytes + return n + + def nbrPatches(self): + n = 0 + for t in self.times(): + for lvl in self.levels(t).values(): + n += len(lvl.patches) + return n + @property def sim(self): if self._sim: diff --git a/pyphare/pyphare/pharesee/run.py b/pyphare/pyphare/pharesee/run.py index ac6d4ac69..0efa03a21 100644 --- a/pyphare/pyphare/pharesee/run.py +++ b/pyphare/pyphare/pharesee/run.py @@ -4,6 +4,7 @@ from .hierarchy.hierarchy_utils import compute_hier_from from .hierarchy.hierarchy_utils import flat_finest_field from .hierarchy import hierarchy_from +from pyphare.core.phare_utilities import listify from pyphare.logger import getLogger @@ -284,12 +285,16 @@ def __init__(self, path, single_hier_for_all_quantities=False): self.single_hier_for_all_quantities = single_hier_for_all_quantities self.hier = None # only used if single_hier_for_all_quantities == True - def _get_hierarchy(self, time, filename, hier=None): - t = "{:.10f}".format(time) + def _get_hierarchy(self, times, filename, hier=None, **kwargs): + times = listify(times) + times = [f"{t:.10f}" for t in times] def _get_hier(h): return hierarchy_from( - h5_filename=os.path.join(self.path, filename), time=t, hier=h + h5_filename=os.path.join(self.path, filename), + times=times, + hier=h, + **kwargs, ) if self.single_hier_for_all_quantities: @@ -297,6 +302,7 @@ def _get_hier(h): return self.hier return _get_hier(hier) + # TODO maybe transform that so multiple times can be accepted def _get(self, hierarchy, time, merged, interp): """ if merged=True, will return an interpolator and a tuple of 1d arrays @@ -322,93 +328,95 @@ def _get(self, hierarchy, time, merged, interp): else: return hierarchy - def GetTags(self, time, merged=False): - hier = self._get_hierarchy(time, "tags.h5") + def GetTags(self, time, merged=False, **kwargs): + hier = self._get_hierarchy(time, "tags.h5", **kwargs) return self._get(hier, time, merged, "nearest") - def GetB(self, time, merged=False, interp="nearest"): - hier = self._get_hierarchy(time, "EM_B.h5") + def GetB(self, time, merged=False, interp="nearest", **kwargs): + hier = self._get_hierarchy(time, "EM_B.h5", **kwargs) return self._get(hier, time, merged, interp) - def GetE(self, time, merged=False, interp="nearest"): - hier = self._get_hierarchy(time, "EM_E.h5") + def GetE(self, time, merged=False, interp="nearest", **kwargs): + hier = self._get_hierarchy(time, "EM_E.h5", **kwargs) return self._get(hier, time, merged, interp) - def GetMassDensity(self, time, merged=False, interp="nearest"): - hier = self._get_hierarchy(time, "ions_mass_density.h5") + def GetMassDensity(self, time, merged=False, interp="nearest", **kwargs): + hier = self._get_hierarchy(time, "ions_mass_density.h5", **kwargs) return self._get(hier, time, merged, interp) - def GetNi(self, time, merged=False, interp="nearest"): - hier = self._get_hierarchy(time, "ions_density.h5") + def GetNi(self, time, merged=False, interp="nearest", **kwargs): + hier = self._get_hierarchy(time, "ions_density.h5", **kwargs) return self._get(hier, time, merged, interp) - def GetN(self, time, pop_name, merged=False, interp="nearest"): - hier = self._get_hierarchy(time, f"ions_pop_{pop_name}_density.h5") + def GetN(self, time, pop_name, merged=False, interp="nearest", **kwargs): + hier = self._get_hierarchy(time, f"ions_pop_{pop_name}_density.h5", **kwargs) return self._get(hier, time, merged, interp) - def GetVi(self, time, merged=False, interp="nearest"): - hier = self._get_hierarchy(time, "ions_bulkVelocity.h5") + def GetVi(self, time, merged=False, interp="nearest", **kwargs): + hier = self._get_hierarchy(time, "ions_bulkVelocity.h5", **kwargs) return self._get(hier, time, merged, interp) - def GetFlux(self, time, pop_name, merged=False, interp="nearest"): - hier = self._get_hierarchy(time, f"ions_pop_{pop_name}_flux.h5") + def GetFlux(self, time, pop_name, merged=False, interp="nearest", **kwargs): + hier = self._get_hierarchy(time, f"ions_pop_{pop_name}_flux.h5", **kwargs) return self._get(hier, time, merged, interp) - def GetPressure(self, time, pop_name, merged=False, interp="nearest"): - M = self._get_hierarchy(time, f"ions_pop_{pop_name}_momentum_tensor.h5") - V = self.GetFlux(time, pop_name) - N = self.GetN(time, pop_name) + def GetPressure(self, time, pop_name, merged=False, interp="nearest", **kwargs): + M = self._get_hierarchy( + time, f"ions_pop_{pop_name}_momentum_tensor.h5", **kwargs + ) + V = self.GetFlux(time, pop_name, **kwargs) + N = self.GetN(time, pop_name, **kwargs) P = compute_hier_from( _compute_pop_pressure, (M, V, N), popname=pop_name, - mass=self.GetMass(pop_name), + mass=self.GetMass(pop_name, **kwargs), ) return self._get(P, time, merged, interp) - def GetPi(self, time, merged=False, interp="nearest"): - M = self._get_hierarchy(time, f"ions_momentum_tensor.h5") - massDensity = self.GetMassDensity(time) - Vi = self._get_hierarchy(time, f"ions_bulkVelocity.h5") + def GetPi(self, time, merged=False, interp="nearest", **kwargs): + M = self._get_hierarchy(time, f"ions_momentum_tensor.h5", **kwargs) + massDensity = self.GetMassDensity(time, **kwargs) + Vi = self._get_hierarchy(time, f"ions_bulkVelocity.h5", **kwargs) Pi = compute_hier_from(_compute_pressure, (M, massDensity, Vi)) return self._get(Pi, time, merged, interp) - def GetJ(self, time, merged=False, interp="nearest"): - B = self.GetB(time) + def GetJ(self, time, merged=False, interp="nearest", **kwargs): + B = self.GetB(time, **kwargs) J = compute_hier_from(_compute_current, B) return self._get(J, time, merged, interp) - def GetDivB(self, time, merged=False, interp="nearest"): - B = self.GetB(time) + def GetDivB(self, time, merged=False, interp="nearest", **kwargs): + B = self.GetB(time, **kwargs) db = compute_hier_from(_compute_divB, B) return self._get(db, time, merged, interp) - def GetRanks(self, time, merged=False, interp="nearest"): + def GetRanks(self, time, merged=False, interp="nearest", **kwargs): """ returns a hierarchy of MPI ranks takes the information from magnetic field diagnostics arbitrarily this fails if the magnetic field is not written and could at some point be replace by a search of any available diag at the requested time. """ - B = self.GetB(time) + B = self.GetB(time, **kwargs) ranks = compute_hier_from(_get_rank, B) return self._get(ranks, time, merged, interp) - def GetParticles(self, time, pop_name, hier=None): + def GetParticles(self, time, pop_name, hier=None, **kwargs): def filename(name): return f"ions_pop_{name}_domain.h5" if isinstance(pop_name, (list, tuple)): for pop in pop_name: - hier = self._get_hierarchy(time, filename(pop), hier=hier) + hier = self._get_hierarchy(time, filename(pop), hier=hier, **kwargs) return hier - return self._get_hierarchy(time, filename(pop_name), hier=hier) + return self._get_hierarchy(time, filename(pop_name), hier=hier, **kwargs) - def GetParticleCount(self, time): - c = self._get_hierarchy(time, "particle_count.h5") + def GetParticleCount(self, time, **kwargs): + c = self._get_hierarchy(time, "particle_count.h5", **kwargs) return c - def GetMass(self, pop_name): + def GetMass(self, pop_name, **kwargs): list_of_qty = ["density", "flux", "domain", "levelGhost", "patchGhost"] list_of_mass = [] @@ -424,7 +432,7 @@ def GetMass(self, pop_name): return list_of_mass[0] - def GetDomainSize(self): + def GetDomainSize(self, **kwargs): h5_filename = "EM_B.h5" # _____ TODO : could be another file import h5py @@ -455,7 +463,7 @@ def GetDl(self, level="finest", time=None): if time is None: time = float(list(data_file[h5_time_grp_key].keys())[0]) - hier = self._get_hierarchy(time, h5_filename) + hier = self._get_hierarchy(time, h5_filename, **kwargs) if level == "finest": level = hier.finest_level(time) @@ -465,7 +473,7 @@ def GetDl(self, level="finest", time=None): return root_cell_width / fac - def GetAllAvailableQties(self, time=0, pops=[]): + def GetAllAvailableQties(self, time=0, pops=[], **kwargs): assert self.single_hier_for_all_quantities # can't work otherwise def _try(fn, *args, **kwargs): @@ -475,14 +483,14 @@ def _try(fn, *args, **kwargs): # normal to not have a diagnostic if not requested logger.debug(f"No file for function {fn.__name__}") - _try(self.GetParticles, time, pops) - _try(self.GetB, time) - _try(self.GetE, time) - _try(self.GetNi, time) - _try(self.GetVi, time) + _try(self.GetParticles, time, pops, **kwargs) + _try(self.GetB, time, **kwargs) + _try(self.GetE, time, **kwargs) + _try(self.GetNi, time, **kwargs) + _try(self.GetVi, time, **kwargs) for pop in pops: - _try(self.GetFlux, time, pop) - _try(self.GetN, time, pop) + _try(self.GetFlux, time, pop, **kwargs) + _try(self.GetN, time, pop, **kwargs) return self.hier From 4c5f5d8e032f2d52b4dbdb2a9c23f91c70a65cd7 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Tue, 4 Jun 2024 14:37:29 +0200 Subject: [PATCH 14/18] up --- pyphare/pyphare/pharesee/hierarchy/fromh5.py | 22 ++++++++++--------- .../pyphare/pharesee/hierarchy/hierarchy.py | 4 +++- tests/functional/tdtagged/td1dtagged.py | 2 +- 3 files changed, 16 insertions(+), 12 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/fromh5.py b/pyphare/pyphare/pharesee/hierarchy/fromh5.py index a8be6fd1f..e8e00353e 100644 --- a/pyphare/pyphare/pharesee/hierarchy/fromh5.py +++ b/pyphare/pyphare/pharesee/hierarchy/fromh5.py @@ -203,13 +203,12 @@ def add_time_from_h5(hier, filepath, time, **kwargs): # we may have a different selection box for that time as for already existing times # but we need to keep them, per time + h5f = h5py.File(filepath, "r") selection_box = kwargs.get("selection_box", None) if hier.has_time(time): - raise ValueError("time already exists in hierarchy") + add_data_from_h5(hier, filepath, time) - patch_levels, h5f = patch_levels_from_h5( - filepath, time, selection_box=selection_box - ) + patch_levels = patch_levels_from_h5(h5f, time, selection_box=selection_box) hier.add_time(time, patch_levels, h5f, selection_box=selection_box) @@ -226,10 +225,14 @@ def add_data_from_h5(hier, filepath, time): if not hier.has_time(time): raise ValueError("time does not exist in hierarchy") + h5f = h5py.File(filepath, "r") + # force using the hierarchy selection box at that time if existing - patch_levels, h5f = patch_levels_from_h5( - filepath, time, selection_box=hier.selection_box[time] - ) + if hier.selection_box is not None: + selection_box = hier.selection_box[time] + else: + selection_box = None + patch_levels = patch_levels_from_h5(h5f, time, selection_box=selection_box) for ilvl, lvl in hier.levels(time).items(): for ip, patch in enumerate(lvl.patches): @@ -285,9 +288,8 @@ def hierarchy_fromh5(h5_filename, time=None, hier=None, silent=True, **kwargs): add_time_from_h5(h, h5_filename, time, **kwargs) if load_one_time(time, hier): - return add_time_from_h5(hier, h5_filename, time, **kwargs) - - return add_data_from_h5(hier, h5_filename, time) + assert len(time) == 1 + return add_time_from_h5(hier, h5_filename, time[0], **kwargs) def hierarchy_fromh5_(h5_filename, time, hier, silent=True): diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index ba3b6fd48..fb63d2a47 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -58,6 +58,8 @@ def __init__( self.data_files.update({data_files.filename: data_files}) else: self.data_files = {data_files.filename: data_files} + else: + self.data_files = {} self.update() @@ -79,7 +81,7 @@ def update(self): new_lvls, self.domain_box, selection_box=self.domain_box, - time=time, + times=time, data_files=self.data_files, ) else: diff --git a/tests/functional/tdtagged/td1dtagged.py b/tests/functional/tdtagged/td1dtagged.py index b418a756c..5b63273ab 100644 --- a/tests/functional/tdtagged/td1dtagged.py +++ b/tests/functional/tdtagged/td1dtagged.py @@ -263,7 +263,7 @@ def get_time(path, time): time = "{:.10f}".format(time) from pyphare.pharesee.hierarchy import hierarchy_from - return hierarchy_from(h5_filename=path + "/ions_pop_protons_domain.h5", time=time) + return hierarchy_from(h5_filename=path + "/ions_pop_protons_domain.h5", times=time) def post_advance(new_time): From fc8d6f9ed94a923a56c244a815d0dfbbaa24d44b Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Tue, 4 Jun 2024 15:57:51 +0200 Subject: [PATCH 15/18] up --- pyphare/pyphare/pharesee/hierarchy/fromh5.py | 176 ++---------------- .../pyphare/pharesee/hierarchy/hierarchy.py | 1 + pyphare/pyphare/pharesee/run.py | 2 +- 3 files changed, 22 insertions(+), 157 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/fromh5.py b/pyphare/pyphare/pharesee/hierarchy/fromh5.py index e8e00353e..0ef71bb7e 100644 --- a/pyphare/pyphare/pharesee/hierarchy/fromh5.py +++ b/pyphare/pyphare/pharesee/hierarchy/fromh5.py @@ -111,10 +111,13 @@ def h5_filename_from(diagInfo): return (diagInfo.quantity + ".h5").replace("/", "_")[1:] -def get_times_from_h5(filepath): +def get_times_from_h5(filepath, as_float=True): f = h5py.File(filepath, "r") - times = np.array(sorted([float(s) for s in list(f["t"].keys())])) + if as_float: + times = np.array(sorted([float(s) for s in list(f["t"].keys())])) + else: + times = list(f["t"].keys()) f.close() return times @@ -208,9 +211,9 @@ def add_time_from_h5(hier, filepath, time, **kwargs): if hier.has_time(time): add_data_from_h5(hier, filepath, time) - patch_levels = patch_levels_from_h5(h5f, time, selection_box=selection_box) - - hier.add_time(time, patch_levels, h5f, selection_box=selection_box) + else: + patch_levels = patch_levels_from_h5(h5f, time, selection_box=selection_box) + hier.add_time(time, patch_levels, h5f, selection_box=selection_box) return hier @@ -251,6 +254,8 @@ def new_from_h5(filepath, times, **kwargs): h5f = h5py.File(filepath, "r") for time in times: + if isinstance(time, float): + time = f"{time:.10f}" patch_levels = patch_levels_from_h5(h5f, time, selection_box=selection_box) patch_levels_per_time.append(patch_levels) @@ -283,161 +288,20 @@ def hierarchy_fromh5(h5_filename, time=None, hier=None, silent=True, **kwargs): if create_from_all_times(time, hier): times = get_times_from_h5(h5_filename) - h = new_from_h5(h5_filename, times[0], **kwargs) - for time in times[1:]: - add_time_from_h5(h, h5_filename, time, **kwargs) + print("creating from all times : ", times, h5_filename) + return new_from_h5(h5_filename, times, **kwargs) if load_one_time(time, hier): assert len(time) == 1 + print("load one time : ", time, hier, h5_filename) return add_time_from_h5(hier, h5_filename, time[0], **kwargs) - -def hierarchy_fromh5_(h5_filename, time, hier, silent=True): - - data_file = h5py.File(h5_filename, "r") - basename = os.path.basename(h5_filename) - - root_cell_width = np.asarray(data_file.attrs["cell_width"]) - interp = data_file.attrs["interpOrder"] - domain_box = Box( - [0] * len(data_file.attrs["domain_box"]), data_file.attrs["domain_box"] - ) - - if create_from_all_times(time, hier): - # first create from first time - # then add all other times - if not silent: - print("creating hierarchy from all times in file") - times = list(data_file[h5_time_grp_key].keys()) - hier = hierarchy_fromh5(h5_filename, time=times[0], hier=hier, silent=silent) - if len(times) > 1: - for t in times[1:]: - hierarchy_fromh5(h5_filename, t, hier, silent=silent) - return hier - - if create_from_one_time(time, hier): - if not silent: - print("creating hierarchy from time {}".format(time)) - t = time - - h5_time_grp = data_file[h5_time_grp_key][time] - patch_levels = {} - - for plvl_key in h5_time_grp.keys(): - h5_patch_lvl_grp = h5_time_grp[plvl_key] - ilvl = int(plvl_key[2:]) - lvl_cell_width = root_cell_width / refinement_ratio**ilvl - patches = {} - - if ilvl not in patches: - patches[ilvl] = [] - - for pkey in h5_patch_lvl_grp.keys(): - h5_patch_grp = h5_patch_lvl_grp[pkey] - - patch_datas = {} - layout = make_layout(h5_patch_grp, lvl_cell_width, interp) - if patch_has_datasets(h5_patch_grp): - add_to_patchdata(patch_datas, h5_patch_grp, basename, layout) - - patches[ilvl].append( - Patch( - patch_datas, - h5_patch_grp.name.split("/")[-1], - layout=layout, - attrs={k: v for k, v in h5_patch_grp.attrs.items()}, - ) - ) - - patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl]) - - diag_hier = PatchHierarchy( - patch_levels, domain_box, refinement_ratio, t, data_files=data_file - ) - - return diag_hier - - if load_one_time(time, hier): - if not silent: - print("loading data at time {} into existing hierarchy".format(time)) - h5_time_grp = data_file[h5_time_grp_key][time] - t = time - - if t in hier.time_hier: - if not silent: - print("time already exist, adding data...") - - # time already exists in the hierarchy - # all we need to do is adding the data - # as patchDatas in the appropriate patches - # and levels, if data compatible with hierarchy - - patch_levels = hier.time_hier[t] - - for plvl_key in h5_time_grp.keys(): - ilvl = int(plvl_key[2:]) - lvl_cell_width = root_cell_width / refinement_ratio**ilvl - - for ipatch, pkey in enumerate(h5_time_grp[plvl_key].keys()): - h5_patch_grp = h5_time_grp[plvl_key][pkey] - - if patch_has_datasets(h5_patch_grp): - hier_patch = patch_levels[ilvl].patches[ipatch] - origin = h5_time_grp[plvl_key][pkey].attrs["origin"] - upper = h5_time_grp[plvl_key][pkey].attrs["upper"] - lower = h5_time_grp[plvl_key][pkey].attrs["lower"] - file_patch_box = Box(lower, upper) - - assert file_patch_box == hier_patch.box - assert (abs(origin - hier_patch.origin) < 1e-6).all() - assert (abs(lvl_cell_width - hier_patch.dl) < 1e-6).all() - - layout = make_layout(h5_patch_grp, lvl_cell_width, interp) - add_to_patchdata( - hier_patch.patch_datas, h5_patch_grp, basename, layout - ) - - return hier - - if not silent: - print("adding data to new time") - # time does not exist in the hierarchy - # we have to create a brand new set of patchLevels - # containing patches, and load data in their patchdatas - - patch_levels = {} - - for plvl_key in h5_time_grp.keys(): - ilvl = int(plvl_key[2:]) - - lvl_cell_width = root_cell_width / refinement_ratio**ilvl - lvl_patches = [] - - for ipatch, pkey in enumerate(h5_time_grp[plvl_key].keys()): - h5_patch_grp = h5_time_grp[plvl_key][pkey] - - layout = make_layout(h5_patch_grp, lvl_cell_width, interp) - patch_datas = {} - if patch_has_datasets(h5_patch_grp): - add_to_patchdata(patch_datas, h5_patch_grp, basename, layout) - lvl_patches.append( - Patch( - patch_datas, - h5_patch_grp.name.split("/")[-1], - layout=layout, - attrs={k: v for k, v in h5_patch_grp.attrs.items()}, - ) - ) - - patch_levels[ilvl] = PatchLevel(ilvl, lvl_patches) - - hier.time_hier[t] = patch_levels - return hier - if load_all_times(time, hier): - if not silent: - print("loading all times in existing hier") - for time in data_file[h5_time_grp_key].keys(): - hier = hierarchy_fromh5(h5_filename, time, hier, silent=silent) - + print("load all times in existing : ", time, h5_filename) + times = get_times_from_h5(h5_filename, as_float=False) + for t in times: + add_time_from_h5(hier, h5_filename, t, **kwargs) return hier + + print("SHOULD NOT BE THERE : ", time, hier) + assert False diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index fb63d2a47..90d2ee864 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -307,6 +307,7 @@ def level_domain_box(self, level_number): def __str__(self): s = "Hierarchy: \n" for t, patch_levels in self.time_hier.items(): + s = s + "Time {}\n".format(t) for ilvl, lvl in patch_levels.items(): s = s + "Level {}\n".format(ilvl) for ip, patch in enumerate(lvl.patches): diff --git a/pyphare/pyphare/pharesee/run.py b/pyphare/pyphare/pharesee/run.py index 0efa03a21..52527b3bd 100644 --- a/pyphare/pyphare/pharesee/run.py +++ b/pyphare/pyphare/pharesee/run.py @@ -463,7 +463,7 @@ def GetDl(self, level="finest", time=None): if time is None: time = float(list(data_file[h5_time_grp_key].keys())[0]) - hier = self._get_hierarchy(time, h5_filename, **kwargs) + hier = self._get_hierarchy(time, h5_filename) if level == "finest": level = hier.finest_level(time) From b89608fb8cd93a516503562b4efa62a0a0d6b557 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Wed, 5 Jun 2024 09:34:45 +0200 Subject: [PATCH 16/18] fix randerror and times arg --- tests/simulator/refinement/test_2d_2_core.py | 6 +++--- tests/simulator/test_advance.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/simulator/refinement/test_2d_2_core.py b/tests/simulator/refinement/test_2d_2_core.py index a50a5c86f..d7ceae427 100644 --- a/tests/simulator/refinement/test_2d_2_core.py +++ b/tests/simulator/refinement/test_2d_2_core.py @@ -143,8 +143,8 @@ def get_time(path, time=None, datahier=None): time = "{:.10f}".format(time) from pyphare.pharesee.hierarchy import hierarchy_from - datahier = hierarchy_from(h5_filename=path + "/EM_E.h5", time=time, hier=datahier) - datahier = hierarchy_from(h5_filename=path + "/EM_B.h5", time=time, hier=datahier) + datahier = hierarchy_from(h5_filename=path + "/EM_E.h5", times=time, hier=datahier) + datahier = hierarchy_from(h5_filename=path + "/EM_B.h5", times=time, hier=datahier) return datahier @@ -188,7 +188,7 @@ def main(): import random startMPI() - rando = random.randint(0, 1e10) + rando = random.randint(0, int(1e10)) refinement_boxes = {"L0": {"B0": [(10, 10), (14, 14)]}} diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py index 4d07aa147..c06f0729f 100644 --- a/tests/simulator/test_advance.py +++ b/tests/simulator/test_advance.py @@ -712,7 +712,7 @@ def _test_field_level_ghosts_via_subcycles_and_coarser_interpolation( import random - rando = random.randint(0, 1e10) + rando = random.randint(0, int(1e10)) def _getHier(diag_dir, boxes=[]): return self.getHierarchy( From d42ca34caadbd5ed1aae812cb82f127232c67e2e Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Wed, 5 Jun 2024 17:19:52 +0200 Subject: [PATCH 17/18] small fixes and add hierarchy pharesee tests --- pyphare/pyphare/pharesee/hierarchy/fromh5.py | 17 +- .../pyphare/pharesee/hierarchy/hierarchy.py | 2 +- pyphare/pyphare/pharesee/run.py | 19 +- .../test_pharesee/CMakeLists.txt | 1 + .../test_pharesee/test_hierarchy.py | 233 ++++++++++++++++++ 5 files changed, 265 insertions(+), 7 deletions(-) create mode 100644 pyphare/pyphare_tests/test_pharesee/test_hierarchy.py diff --git a/pyphare/pyphare/pharesee/hierarchy/fromh5.py b/pyphare/pyphare/pharesee/hierarchy/fromh5.py index 0ef71bb7e..d84b9d981 100644 --- a/pyphare/pyphare/pharesee/hierarchy/fromh5.py +++ b/pyphare/pyphare/pharesee/hierarchy/fromh5.py @@ -7,7 +7,12 @@ from ..particles import Particles from .hierarchy import PatchHierarchy from ...core.box import Box -from ...core.phare_utilities import refinement_ratio +from ...core.phare_utilities import ( + refinement_ratio, + none_iterable, + all_iterables, + listify, +) from ...core.gridlayout import GridLayout from .hierarchy_utils import field_qties import h5py @@ -248,15 +253,19 @@ def new_from_h5(filepath, times, **kwargs): # create a patchhierarchy from a given time and optional selection box # loads all datasets from the filepath h5 file as patchdatas - selection_box = kwargs.get("selection_box", None) + # we authorize user to pass only one selection box for all times + # but in this case they're all the same + selection_box = kwargs.get("selection_box", [None] * len(times)) + if none_iterable(selection_box) and all_iterables(times): + selection_box = [selection_box] * len(times) patch_levels_per_time = [] h5f = h5py.File(filepath, "r") - for time in times: + for it, time in enumerate(times): if isinstance(time, float): time = f"{time:.10f}" - patch_levels = patch_levels_from_h5(h5f, time, selection_box=selection_box) + patch_levels = patch_levels_from_h5(h5f, time, selection_box=selection_box[it]) patch_levels_per_time.append(patch_levels) dim = len(h5f.attrs["domain_box"]) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index 90d2ee864..f263a14d5 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -330,7 +330,7 @@ def has_file(self, filename): return filename in self.data_files def times(self): - return np.sort(np.asarray(list(self.time_hier.keys()))) + return np.sort(np.asarray(list(self.time_hier.keys()), dtype=np.float32)) def plot_patches(self, save=False): fig, ax = plt.subplots(figsize=(10, 3)) diff --git a/pyphare/pyphare/pharesee/run.py b/pyphare/pyphare/pharesee/run.py index 52527b3bd..3bbd14738 100644 --- a/pyphare/pyphare/pharesee/run.py +++ b/pyphare/pyphare/pharesee/run.py @@ -280,14 +280,27 @@ def make_interpolator(data, coords, interp, domain, dl, qty, nbrGhosts): class Run: + def __init__(self, path, single_hier_for_all_quantities=False): + import glob + self.path = path self.single_hier_for_all_quantities = single_hier_for_all_quantities self.hier = None # only used if single_hier_for_all_quantities == True + self.available_diags = glob.glob(os.path.join(self.path, "*.h5")) + if len(self.available_diags) == 0: + raise FileNotFoundError(f"No diagnostic files found in {self.path}") def _get_hierarchy(self, times, filename, hier=None, **kwargs): + from pyphare.core.box import Box + times = listify(times) times = [f"{t:.10f}" for t in times] + if "selection_box" in kwargs: + if isinstance(kwargs["selection_box"], tuple): + lower = kwargs["selection_box"][:2] + upper = kwargs["selection_box"][2:] + kwargs["selection_box"] = Box(lower, upper) def _get_hier(h): return hierarchy_from( @@ -452,10 +465,12 @@ def GetDl(self, level="finest", time=None): :param level: the level at which get the associated grid size :param time: the time because level depends on it """ + import glob h5_time_grp_key = "t" - h5_filename = "EM_B.h5" # _____ TODO : could be another file - + files = glob.glob(os.path.join(self.path, "*.h5")) + any_file = files[0] + h5_filename = any_file import h5py data_file = h5py.File(os.path.join(self.path, h5_filename), "r") diff --git a/pyphare/pyphare_tests/test_pharesee/CMakeLists.txt b/pyphare/pyphare_tests/test_pharesee/CMakeLists.txt index 133d31b65..b40969d48 100644 --- a/pyphare/pyphare_tests/test_pharesee/CMakeLists.txt +++ b/pyphare/pyphare_tests/test_pharesee/CMakeLists.txt @@ -5,5 +5,6 @@ project(test-pharesee) add_python3_test(test-pharesee-geometry_1d test_geometry.py ${PROJECT_SOURCE_DIR}) add_python3_test(test-pharesee-geometry_2d test_geometry_2d.py ${PROJECT_SOURCE_DIR}) +phare_mpi_python3_exec(11 2 test-pharesee-hierarchy test_hierarchy.py ${PROJECT_SOURCE_DIR}) diff --git a/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py new file mode 100644 index 000000000..65bf02a42 --- /dev/null +++ b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py @@ -0,0 +1,233 @@ +import unittest +from ddt import ddt, data, unpack +import numpy as np + +from pyphare.simulator.simulator import Simulator +from pyphare.pharesee.run import Run +from pyphare.pharesee.hierarchy import PatchHierarchy +from pyphare.core.box import Box + +diag_outputs = "phare_outputs/" +time_step_nbr = 6000 +time_step = 0.005 +final_time = time_step * time_step_nbr +dt = 10 * time_step +nt = final_time / dt + 1 +timestamps = dt * np.arange(nt) + + +@ddt +class PatchHierarchyTest(unittest.TestCase): + + @classmethod + def setUpClass(cls): + from pyphare.pharein import global_vars + import pyphare.pharein as ph + + global_vars.sim = None + + def config(): + sim = ph.Simulation( + time_step_nbr=time_step_nbr, + time_step=time_step, + cells=(86, 86), + dl=(0.3, 0.3), + refinement="tagging", + max_nbr_levels=2, + hyper_resistivity=0.005, + resistivity=0.001, + diag_options={ + "format": "phareh5", + "options": {"dir": diag_outputs, "mode": "overwrite"}, + }, + ) + + def density(x, y): + L = sim.simulation_domain()[1] + return ( + 0.2 + + 1.0 / np.cosh((y - L * 0.3) / 0.5) ** 2 + + 1.0 / np.cosh((y - L * 0.7) / 0.5) ** 2 + ) + + def S(y, y0, l): + return 0.5 * (1.0 + np.tanh((y - y0) / l)) + + def by(x, y): + Lx = sim.simulation_domain()[0] + Ly = sim.simulation_domain()[1] + w1 = 0.2 + w2 = 1.0 + x0 = x - 0.5 * Lx + y1 = y - 0.3 * Ly + y2 = y - 0.7 * Ly + w3 = np.exp(-(x0 * x0 + y1 * y1) / (w2 * w2)) + w4 = np.exp(-(x0 * x0 + y2 * y2) / (w2 * w2)) + w5 = 2.0 * w1 / w2 + return (w5 * x0 * w3) + (-w5 * x0 * w4) + + def bx(x, y): + Lx = sim.simulation_domain()[0] + Ly = sim.simulation_domain()[1] + w1 = 0.2 + w2 = 1.0 + x0 = x - 0.5 * Lx + y1 = y - 0.3 * Ly + y2 = y - 0.7 * Ly + w3 = np.exp(-(x0 * x0 + y1 * y1) / (w2 * w2)) + w4 = np.exp(-(x0 * x0 + y2 * y2) / (w2 * w2)) + w5 = 2.0 * w1 / w2 + v1 = -1 + v2 = 1.0 + return ( + v1 + + (v2 - v1) * (S(y, Ly * 0.3, 0.5) - S(y, Ly * 0.7, 0.5)) + + (-w5 * y1 * w3) + + (+w5 * y2 * w4) + ) + + def bz(x, y): + return 0.0 + + def b2(x, y): + return bx(x, y) ** 2 + by(x, y) ** 2 + bz(x, y) ** 2 + + def T(x, y): + K = 1 + temp = 1.0 / density(x, y) * (K - b2(x, y) * 0.5) + assert np.all(temp > 0) + return temp + + def vx(x, y): + return 0.0 + + def vy(x, y): + return 0.0 + + def vz(x, y): + return 0.0 + + def vthx(x, y): + return np.sqrt(T(x, y)) + + def vthy(x, y): + return np.sqrt(T(x, y)) + + def vthz(x, y): + return np.sqrt(T(x, y)) + + vvv = { + "vbulkx": vx, + "vbulky": vy, + "vbulkz": vz, + "vthx": vthx, + "vthy": vthy, + "vthz": vthz, + "nbr_part_per_cell": 100, + } + + ph.MaxwellianFluidModel( + bx=bx, + by=by, + bz=bz, + protons={ + "charge": 1, + "density": density, + **vvv, + "init": {"seed": 12334}, + }, + ) + + ph.ElectronModel(closure="isothermal", Te=0.1) + + for quantity in ["E", "B"]: + ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) + + for quantity in ["density", "bulkVelocity"]: + ph.FluidDiagnostics( + quantity=quantity, + write_timestamps=timestamps, + compute_timestamps=timestamps, + ) + + return sim + + Simulator(config()).run() + + def test_data_is_a_hierarchy(self): + r = Run(diag_outputs) + B = r.GetB(0.0) + self.assertTrue(isinstance(B, PatchHierarchy)) + + def test_can_read_multiple_times(self): + r = Run(diag_outputs) + times = (0.0, 0.1) + B = r.GetB(times) + E = r.GetE(times) + Ni = r.GetNi(times) + Vi = r.GetVi(times) + for hier in (B, E, Ni, Vi): + self.assertEqual(len(hier.times()), 2) + self.assertTrue(np.allclose(hier.times(), times)) + + def test_hierarchy_is_refined(self): + r = Run(diag_outputs) + time = 0.0 + B = r.GetB(time) + self.assertEqual(len(B.levels()), self.levelNbrs()) + self.assertEqual(len(B.levels()), 2) + self.assertEqual(len(B.levels(time)), 2) + self.assertEqual(len(B.levels(time)), self.levelNbrs(time)) + + def test_can_get_nbytes(self): + r = Run(diag_outputs) + time = 0.0 + B = r.GetB(time) + self.assertTrue(B.nbytes() > 0) + + def test_hierarchy_has_patches(self): + r = Run(diag_outputs) + time = 0.0 + B = r.GetB(time) + self.assertTrue(B.nbrPatches() > 0) + + def test_access_patchdatas_as_hierarchies(self): + r = Run(diag_outputs) + time = 0.0 + B = r.GetB(time) + self.assertTrue(isinstance(B.Bx, PatchHierarchy)) + self.assertTrue(isinstance(B.By, PatchHierarchy)) + self.assertTrue(isinstance(B.Bz, PatchHierarchy)) + + def test_partial_domain_hierarchies(self): + import matplotlib.pyplot as plt + from matplotlib.patches import Rectangle + + r = Run(diag_outputs) + time = 0.0 + box = Box((10, 5), (18, 12.5)) + B = r.GetB(time) + Bpartial = r.GetB(time, selection_box=box) + + fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5)) + B.Bx.plot(plot_patches=True, ax=ax1) + Bpartial.Bx.plot(plot_patches=True, ax=ax2) + ax2.add_patch( + Rectangle( + box.lower, + box.shape[0], + box.shape[1], + ec="w", + fc="none", + lw=3, + ) + ) + fig.set_title(f"{self.id()}") + fig.savefig(f"{self.id()}.png") + + self.assertTrue(isinstance(Bpartial, PatchHierarchy)) + self.assertTrue(Bpartial.nbrPatches() < B.nbrPatches()) + + +if __name__ == "__main__": + unittest.main() From ad2781ce97eaa3783f2128a78a13579cb88b5d35 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Wed, 5 Jun 2024 18:55:13 +0200 Subject: [PATCH 18/18] up --- pyphare/pyphare/pharesee/hierarchy/hierarchy.py | 3 ++- pyphare/pyphare/pharesee/run.py | 2 -- pyphare/pyphare_tests/test_pharesee/test_hierarchy.py | 2 +- tests/simulator/refinement/test_2d_10_core.py | 4 ++-- 4 files changed, 5 insertions(+), 6 deletions(-) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index f263a14d5..844b1e5e6 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -330,7 +330,8 @@ def has_file(self, filename): return filename in self.data_files def times(self): - return np.sort(np.asarray(list(self.time_hier.keys()), dtype=np.float32)) + # return np.sort(np.asarray(list(self.time_hier.keys()), dtype=np.float32)) + return np.sort(np.asarray(list(self.time_hier.keys()))) def plot_patches(self, save=False): fig, ax = plt.subplots(figsize=(10, 3)) diff --git a/pyphare/pyphare/pharesee/run.py b/pyphare/pyphare/pharesee/run.py index 3bbd14738..baf5e3a2e 100644 --- a/pyphare/pyphare/pharesee/run.py +++ b/pyphare/pyphare/pharesee/run.py @@ -288,8 +288,6 @@ def __init__(self, path, single_hier_for_all_quantities=False): self.single_hier_for_all_quantities = single_hier_for_all_quantities self.hier = None # only used if single_hier_for_all_quantities == True self.available_diags = glob.glob(os.path.join(self.path, "*.h5")) - if len(self.available_diags) == 0: - raise FileNotFoundError(f"No diagnostic files found in {self.path}") def _get_hierarchy(self, times, filename, hier=None, **kwargs): from pyphare.core.box import Box diff --git a/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py index 65bf02a42..0eda79638 100644 --- a/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py +++ b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py @@ -168,7 +168,7 @@ def test_can_read_multiple_times(self): Vi = r.GetVi(times) for hier in (B, E, Ni, Vi): self.assertEqual(len(hier.times()), 2) - self.assertTrue(np.allclose(hier.times(), times)) + self.assertTrue(np.allclose(hier.times().astype(np.float32), times)) def test_hierarchy_is_refined(self): r = Run(diag_outputs) diff --git a/tests/simulator/refinement/test_2d_10_core.py b/tests/simulator/refinement/test_2d_10_core.py index 2913bf569..c548db4ff 100644 --- a/tests/simulator/refinement/test_2d_10_core.py +++ b/tests/simulator/refinement/test_2d_10_core.py @@ -120,8 +120,8 @@ def get_time(path, time=None, datahier=None): time = "{:.10f}".format(time) from pyphare.pharesee.hierarchy import hierarchy_from - datahier = hierarchy_from(h5_filename=path + "/EM_E.h5", time=time, hier=datahier) - datahier = hierarchy_from(h5_filename=path + "/EM_B.h5", time=time, hier=datahier) + datahier = hierarchy_from(h5_filename=path + "/EM_E.h5", times=time, hier=datahier) + datahier = hierarchy_from(h5_filename=path + "/EM_B.h5", times=time, hier=datahier) return datahier