diff --git a/pyphare/pyphare/core/box.py b/pyphare/pyphare/core/box.py index fad467dbe..9b59b3a1c 100644 --- a/pyphare/pyphare/core/box.py +++ b/pyphare/pyphare/core/box.py @@ -70,26 +70,26 @@ def copy(self): class nDBox(Box): - def __init__(self, dim, l, u): + def __init__(self, dim, lower, upper): def _get(self, p): return np.asarray([p] * dim) - super().__init__(_get(dim, l), _get(dim, u)) + super().__init__(_get(dim, lower), _get(dim, upper)) class Box1D(nDBox): - def __init__(self, l, u): - super().__init__(1, l, u) + def __init__(self, lower, upper): + super().__init__(1, lower, upper) class Box2D(nDBox): - def __init__(self, l, u): - super().__init__(2, l, u) + def __init__(self, lower, upper): + super().__init__(2, lower, upper) class Box3D(nDBox): - def __init__(self, l, u): - super().__init__(3, l, u) + def __init__(self, lower, upper): + super().__init__(3, lower, upper) def refine(box, ratio): diff --git a/pyphare/pyphare/core/operators.py b/pyphare/pyphare/core/operators.py new file mode 100644 index 000000000..1bd81fb12 --- /dev/null +++ b/pyphare/pyphare/core/operators.py @@ -0,0 +1,141 @@ +import numpy as np + +from pyphare.pharesee.hierarchy import ScalarField, VectorField +from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from +from pyphare.pharesee.hierarchy.hierarchy_utils import rename + + +def _compute_dot_product(patch_datas, **kwargs): + ref_name = next(iter(patch_datas.keys())) + + dset = ( + patch_datas["left_x"][:] * patch_datas["right_x"][:] + + patch_datas["left_y"][:] * patch_datas["right_y"][:] + + patch_datas["left_z"][:] * patch_datas["right_z"][:] + ) + + return ( + {"name": "value", "data": dset, "centering": patch_datas[ref_name].centerings}, + ) + + +def _compute_sqrt(patch_datas, **kwargs): + ref_name = next(iter(patch_datas.keys())) + + dset = np.sqrt(patch_datas["value"][:]) + + return ( + {"name": "value", "data": dset, "centering": patch_datas[ref_name].centerings}, + ) + + +def _compute_cross_product(patch_datas, **kwargs): + ref_name = next(iter(patch_datas.keys())) + + dset_x = ( + patch_datas["left_y"][:] * patch_datas["right_z"][:] + - patch_datas["left_z"][:] * patch_datas["right_y"][:] + ) + dset_y = ( + patch_datas["left_z"][:] * patch_datas["right_x"][:] + - patch_datas["left_x"][:] * patch_datas["right_z"][:] + ) + dset_z = ( + patch_datas["left_x"][:] * patch_datas["right_y"][:] + - patch_datas["left_y"][:] * patch_datas["right_x"][:] + ) + + return ( + {"name": "x", "data": dset_x, "centering": patch_datas[ref_name].centerings}, + {"name": "y", "data": dset_y, "centering": patch_datas[ref_name].centerings}, + {"name": "z", "data": dset_z, "centering": patch_datas[ref_name].centerings}, + ) + + +def _compute_grad(patch_data, **kwargs): + ndim = patch_data["value"].box.ndim + nb_ghosts = kwargs["nb_ghosts"] + ds = patch_data["value"].dataset + + ds_shape = list(ds.shape) + + ds_x = np.full(ds_shape, np.nan) + ds_y = np.full(ds_shape, np.nan) + ds_z = np.full(ds_shape, np.nan) + + grad_ds = np.gradient(ds) + select = tuple([slice(nb_ghosts, -nb_ghosts) for _ in range(ndim)]) + if ndim == 2: + ds_x[select] = np.asarray(grad_ds[0][select]) + ds_y[select] = np.asarray(grad_ds[1][select]) + ds_z[select].fill(0.0) # TODO at 2D, gradient is null in z dir + + else: + raise RuntimeError("dimension not yet implemented") + + return ( + {"name": "x", "data": ds_x, "centering": patch_data["value"].centerings}, + {"name": "y", "data": ds_y, "centering": patch_data["value"].centerings}, + {"name": "z", "data": ds_z, "centering": patch_data["value"].centerings}, + ) + + +def dot(hier_left, hier_right, **kwargs): + if isinstance(hier_left, VectorField) and isinstance(hier_right, VectorField): + names_left = ["left_x", "left_y", "left_z"] + names_right = ["right_x", "right_y", "right_z"] + + else: + raise RuntimeError("type of hierarchy not yet considered") + + hl = rename(hier_left, names_left) + hr = rename(hier_right, names_right) + + h = compute_hier_from( + _compute_dot_product, + (hl, hr), + ) + + return ScalarField(h) + + +def cross(hier_left, hier_right, **kwargs): + if isinstance(hier_left, VectorField) and isinstance(hier_right, VectorField): + names_left = ["left_x", "left_y", "left_z"] + names_right = ["right_x", "right_y", "right_z"] + + else: + raise RuntimeError("type of hierarchy not yet considered") + + hl = rename(hier_left, names_left) + hr = rename(hier_right, names_right) + + h = compute_hier_from( + _compute_cross_product, + (hl, hr), + ) + + return VectorField(h) + + +def sqrt(hier, **kwargs): + h = compute_hier_from( + _compute_sqrt, + hier, + ) + + return ScalarField(h) + + +def modulus(hier): + assert isinstance(hier, VectorField) + + return sqrt(dot(hier, hier)) + + +def grad(hier, **kwargs): + assert isinstance(hier, ScalarField) + nb_ghosts = list(hier.level(0).patches[0].patch_datas.values())[0].ghosts_nbr[0] + h = compute_hier_from(_compute_grad, hier, nb_ghosts=nb_ghosts) + + return VectorField(h) diff --git a/pyphare/pyphare/core/phare_utilities.py b/pyphare/pyphare/core/phare_utilities.py index 27a092129..c6937cc76 100644 --- a/pyphare/pyphare/core/phare_utilities.py +++ b/pyphare/pyphare/core/phare_utilities.py @@ -155,7 +155,8 @@ def run_cli_cmd(cmd, shell=True, capture_output=True, check=False, print_cmd=Fal def print_trace(): - import sys, traceback + import sys + import traceback _, _, tb = sys.exc_info() traceback.print_tb(tb) diff --git a/pyphare/pyphare/cpp/__init__.py b/pyphare/pyphare/cpp/__init__.py index 79d0c8cf5..c66b11d38 100644 --- a/pyphare/pyphare/cpp/__init__.py +++ b/pyphare/pyphare/cpp/__init__.py @@ -22,7 +22,7 @@ def cpp_lib(override=None): return importlib.import_module("pybindlibs.cpp") try: return importlib.import_module("pybindlibs.cpp_dbg") - except ImportError as err: + except ImportError: return importlib.import_module("pybindlibs.cpp") diff --git a/pyphare/pyphare/pharein/__init__.py b/pyphare/pyphare/pharein/__init__.py index 21beb2989..202b13eee 100644 --- a/pyphare/pyphare/pharein/__init__.py +++ b/pyphare/pyphare/pharein/__init__.py @@ -4,6 +4,34 @@ import numpy as np from pyphare.core.phare_utilities import is_scalar +from .uniform_model import UniformModel +from .maxwellian_fluid_model import MaxwellianFluidModel +from .electron_model import ElectronModel +from .diagnostics import ( + FluidDiagnostics, + ElectromagDiagnostics, + ParticleDiagnostics, + MetaDiagnostics, + InfoDiagnostics, +) +from .simulation import ( + Simulation, + serialize as serialize_sim, + deserialize as deserialize_sim, +) +from .load_balancer import LoadBalancer + +__all__ = [ + "UniformModel", + "MaxwellianFluidModel", + "ElectronModel", + "FluidDiagnostics", + "ElectromagDiagnostics", + "ParticleDiagnostics", + "MetaDiagnostics", + "InfoDiagnostics", + "Simulation", +] # This exists to allow a condition variable for when we are running PHARE from C++ via phare-exe # It is configured to "True" in pyphare/pyphare/pharein/init.py::get_user_inputs(jobname) @@ -29,24 +57,6 @@ sys.path = sys.path + pythonpath -from .uniform_model import UniformModel -from .maxwellian_fluid_model import MaxwellianFluidModel -from .electron_model import ElectronModel -from .diagnostics import ( - FluidDiagnostics, - ElectromagDiagnostics, - ParticleDiagnostics, - MetaDiagnostics, - InfoDiagnostics, -) -from .simulation import ( - Simulation, - serialize as serialize_sim, - deserialize as deserialize_sim, -) -from .load_balancer import LoadBalancer - - def NO_GUI(): """prevents issues when command line only and no desktop etc""" import matplotlib as mpl diff --git a/pyphare/pyphare/pharein/diagnostics.py b/pyphare/pyphare/pharein/diagnostics.py index e2f136481..9eb5ce4e3 100644 --- a/pyphare/pyphare/pharein/diagnostics.py +++ b/pyphare/pyphare/pharein/diagnostics.py @@ -147,7 +147,7 @@ def __init__(self, name, **kwargs): self.__extent = None - # if a diag already is registered we just contactenate the timestamps + # if a diag already is registered we just concatenate the timestamps addIt = True registered_diags = global_vars.sim.diagnostics for diagname, diag in registered_diags.items(): @@ -313,7 +313,6 @@ def __init__(self, **kwargs): ) def _setSubTypeAttributes(self, **kwargs): - # domain is good default for users who should not worry about what that means # even less about ghosts... kwargs["quantity"] = kwargs.get("quantity", "domain") diff --git a/pyphare/pyphare/pharein/examples/job.py b/pyphare/pyphare/pharein/examples/job.py index 5b7cb08a2..24ed09f98 100644 --- a/pyphare/pyphare/pharein/examples/job.py +++ b/pyphare/pyphare/pharein/examples/job.py @@ -1,25 +1,20 @@ #!/usr/bin/env python - -from pyphare.pharein import Simulation -from pyphare.pharein import MaxwellianFluidModel -from pyphare.pharein import ElectronModel -from pyphare.pharein import ElectromagDiagnostics -from pyphare.pharein import FluidDiagnostics -from pyphare.pharein import getSimulation +import numpy as np +import pyphare.pharein as ph # ------------------------------------ # configure the simulation # ------------------------------------ -Simulation( +ph.Simulation( time_step_nbr=1000, # number of time steps (not specified if time_step and final_time provided) final_time=1.0, # simulation final time (not specified if time_step and time_step_nbr given) boundary_types="periodic", # boundary condition, string or tuple, length == len(cell) == len(dl) cells=80, # integer or tuple length == dimension dl=0.1, # mesh size of the root level, float or tuple - path="test5" # directory where INI file and diagnostics directories will be + path="test5", # directory where INI file and diagnostics directories will be # time_step = 0.005, # simulation time step (not specified if time_step_nbr and final_time given) # domain_size = 8., # float or tuple, not specified if dl and cells are # interp_order = 1, # interpolation order, [default = 1] can be 1, 2, 3 or 4 @@ -36,7 +31,6 @@ # in the following we use the MaxwellianFluidModel -import numpy as np Te = 0.12 @@ -46,17 +40,17 @@ def n(x): def bx(x): - xmax = getSimulation().simulation_domain()[0] + xmax = ph.getSimulation().simulation_domain()[0] return np.cos(2 * np.pi / xmax * x) -MaxwellianFluidModel(bx=bx, protons={"density": n}, background={}) +ph.MaxwellianFluidModel(bx=bx, protons={"density": n}, background={}) -ElectronModel(closure="isothermal", Te=Te) +ph.ElectronModel(closure="isothermal", Te=Te) -ElectromagDiagnostics( +ph.ElectromagDiagnostics( diag_type="E", # available : ("E", "B") write_every=10, compute_every=5, @@ -66,18 +60,18 @@ def bx(x): ) -FluidDiagnostics( +ph.FluidDiagnostics( diag_type="density", # choose in (rho_s, flux_s) write_every=10, # write on disk every x iterations compute_every=5, # compute diagnostics every x iterations ( x <= write_every) start_iteration=0, # iteration at which diag is enabled last_iteration=990, # iteration at which diag is turned off - population_name="protons" # name of the population for which the diagnostics is made + population_name="protons", # name of the population for which the diagnostics is made # ,path = 'FluidDiagnostics1' # where output files will be written, [default: name] ) -FluidDiagnostics( +ph.FluidDiagnostics( diag_type="bulkVelocity", write_every=10, compute_every=5, @@ -86,7 +80,7 @@ def bx(x): population_name="background", ) -FluidDiagnostics( +ph.FluidDiagnostics( diag_type="density", write_every=10, compute_every=5, @@ -95,7 +89,7 @@ def bx(x): population_name="all", ) -FluidDiagnostics( +ph.FluidDiagnostics( diag_type="flux", write_every=10, compute_every=5, @@ -104,9 +98,13 @@ def bx(x): population_name="background", ) -ElectromagDiagnostics( - diag_type="B", write_every=10, compute_every=5, start_teration=0, last_iteration=990 +ph.ElectromagDiagnostics( + diag_type="B", + write_every=10, + compute_every=5, + start_iteration=0, + last_iteration=990, ) -for item in getSimulation().electrons.dict_path(): +for item in ph.getSimulation().electrons.dict_path(): print(item[0], item[1]) diff --git a/pyphare/pyphare/pharein/maxwellian_fluid_model.py b/pyphare/pyphare/pharein/maxwellian_fluid_model.py index 62a37381c..4e3f678dd 100644 --- a/pyphare/pyphare/pharein/maxwellian_fluid_model.py +++ b/pyphare/pyphare/pharein/maxwellian_fluid_model.py @@ -1,7 +1,7 @@ import numpy as np from pyphare.core import phare_utilities from pyphare.core.box import Box -from pyphare.core.gridlayout import GridLayout, yee_element_is_primal +from pyphare.core.gridlayout import GridLayout from pyphare.pharein import global_vars diff --git a/pyphare/pyphare/pharein/restarts.py b/pyphare/pyphare/pharein/restarts.py index 4d663e8a4..fefab3133 100644 --- a/pyphare/pyphare/pharein/restarts.py +++ b/pyphare/pyphare/pharein/restarts.py @@ -21,14 +21,14 @@ def validate(sim): if not np.all(np.diff(restart_options["elapsed_timestamps"]) >= 0): raise RuntimeError( - f"Error: restart_options elapsed_timestamps not in ascending order)" + "Error: restart_options elapsed_timestamps not in ascending order)" ) # seconds_in_an_hour = 60 ** 2 # for cmp_idx, ref_ts in enumerate(restart_options["elapsed_timestamps"][1:]): # cmp_ts = restart_options["elapsed_timestamps"][cmp_idx] # if ref_ts - cmp_ts < seconds_in_an_hour: - # raise RuntimeError(f"Error: time betweeen restart_options elapsed_timestamps must be at least one hour)") + # raise RuntimeError("Error: time betweeen restart_options elapsed_timestamps must be at least one hour)") if "timestamps" in restart_options: restart_options["timestamps"] = phare_utilities.np_array_ify( @@ -47,7 +47,7 @@ def validate(sim): ) if not np.all(np.diff(timestamps) >= 0): raise RuntimeError( - f"Error: restart_options timestamps not in ascending order)" + "Error: restart_options timestamps not in ascending order)" ) if not np.all( np.abs( @@ -55,7 +55,7 @@ def validate(sim): ) ): raise RuntimeError( - f"Error: restart_options timestamps is inconsistent with simulation.time_step" + "Error: restart_options timestamps is inconsistent with simulation.time_step" ) sim.restart_options["timestamps"] = conserve_existing(sim) diff --git a/pyphare/pyphare/pharein/simulation.py b/pyphare/pyphare/pharein/simulation.py index e66558ac9..b8f0fa787 100644 --- a/pyphare/pyphare/pharein/simulation.py +++ b/pyphare/pyphare/pharein/simulation.py @@ -332,7 +332,7 @@ def check_refinement_boxes(ndim, **kwargs): for box_idx, ref_box in enumerate(boxes): for cmp_box in boxes[box_idx + 1 :]: - if ref_box * cmp_box != None: + if ref_box * cmp_box is not None: raise ValueError( f"Error: Box({ref_box}) overlaps with Box({cmp_box})" ) @@ -352,8 +352,8 @@ def check_refinement_boxes(ndim, **kwargs): if box.ndim != ndim: raise ValueError(f"Box({box}) has incorrect dimensions for simulation") - for l in boxm.refine(box, refinement_ratio).shape: - if (l < smallest_patch_size).any(): + for length in boxm.refine(box, refinement_ratio).shape: + if (length < smallest_patch_size).any(): raise ValueError( "Invalid box incompatible with smallest_patch_size" ) @@ -558,7 +558,7 @@ def check_optional_keywords(**kwargs): def check_resistivity(**kwargs): resistivity = kwargs.get("resistivity", 0.0) if resistivity < 0.0: - raise ValueError(f"Error: resistivity should not be negative") + raise ValueError("Error: resistivity should not be negative") return resistivity @@ -566,7 +566,7 @@ def check_resistivity(**kwargs): def check_hyper_resistivity(**kwargs): hyper_resistivity = kwargs.get("hyper_resistivity", 0.0001) if hyper_resistivity < 0.0: - raise ValueError(f"Error: hyper_resistivity should not be negative") + raise ValueError("Error: hyper_resistivity should not be negative") return hyper_resistivity @@ -675,7 +675,7 @@ def wrapper(simulation_object, **kwargs): ) = check_refinement_boxes(ndim, **kwargs) else: kwargs["max_nbr_levels"] = kwargs.get("max_nbr_levels", None) - assert kwargs["max_nbr_levels"] != None # this needs setting otherwise + assert kwargs["max_nbr_levels"] is not None # this needs setting otherwise kwargs["refinement_boxes"] = None kwargs["resistivity"] = check_resistivity(**kwargs) @@ -936,12 +936,14 @@ def set_electrons(self, electrons): def serialize(sim): # pickle cannot handle simulation objects - import dill, codecs + import dill + import codecs return codecs.encode(dill.dumps(sim), "hex") def deserialize(hex): - import dill, codecs + import dill + import codecs return dill.loads(codecs.decode(hex, "hex")) 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 d3a7fef1a..000000000 --- a/pyphare/pyphare/pharesee/hierarchy.py +++ /dev/null @@ -1,1861 +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 - -from dataclasses import dataclass - - -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 __eq__(self, that): - return self.field_name == that.field_name and self.dataset[:] == that.dataset[:] - - 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() - - def __eq__(self, that): - return self.name == that.name and self.dataset == that.dataset - - -class Patch: - """ - A patch represents a hyper-rectangular region of space - """ - - def __init__(self, patch_datas, patch_id=""): - """ - :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 - - def __str__(self): - return f"Patch: box( {self.box}), id({self.id})" - - def __repr__(self): - return self.__str__() - - def __getitem__(self, key): - return self.patch_datas[key] - - 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) - if len(times) > 1: - for t in times[1:]: - hierarchy_fromh5(h5_filename, t, hier) - 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] - - 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]) - ) - - 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] - - 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)) - - 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) - - 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]] - return hierarchy - - -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 - - -@dataclass -class EqualityReport: - ok: bool - reason: str - - def __bool__(self): - return self.ok - - -def hierarchy_compare(this, that): - if not isinstance(this, PatchHierarchy) or not isinstance(that, PatchHierarchy): - return EqualityReport(False, "class type mismatch") - - if this.ndim != that.ndim or this.domain_box != that.domain_box: - return EqualityReport(False, "dimensional mismatch") - - if this.time_hier.keys() != that.time_hier.keys(): - return EqualityReport(False, "timesteps mismatch") - - for tidx in this.times(): - patch_levels_ref = this.time_hier[tidx] - patch_levels_cmp = that.time_hier[tidx] - - if patch_levels_ref.keys() != patch_levels_cmp.keys(): - return EqualityReport(False, "levels mismatch") - - for level_idx in patch_levels_cmp.keys(): - patch_level_ref = patch_levels_ref[level_idx] - patch_level_cmp = patch_levels_cmp[level_idx] - - for patch_idx in range(len(patch_level_cmp.patches)): - patch_ref = patch_level_ref.patches[patch_idx] - patch_cmp = patch_level_cmp.patches[patch_idx] - - if patch_ref.patch_datas.keys() != patch_cmp.patch_datas.keys(): - print(list(patch_ref.patch_datas.keys())) - print(list(patch_cmp.patch_datas.keys())) - return EqualityReport(False, "data keys mismatch") - - for patch_data_key in patch_ref.patch_datas.keys(): - patch_data_ref = patch_ref.patch_datas[patch_data_key] - patch_data_cmp = patch_cmp.patch_datas[patch_data_key] - - if patch_data_cmp != patch_data_ref: - return EqualityReport( - False, - "data mismatch: " - + type(patch_data_cmp).__name__ - + " " - + type(patch_data_ref).__name__, - ) - - return EqualityReport(True, "OK") diff --git a/pyphare/pyphare/pharesee/hierarchy/__init__.py b/pyphare/pyphare/pharesee/hierarchy/__init__.py new file mode 100644 index 000000000..6bc7e43c8 --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/__init__.py @@ -0,0 +1,42 @@ +from .scalarfield import ScalarField +from .vectorfield import VectorField +from .hierarchy import PatchHierarchy +from pyphare.core.phare_utilities import listify + +__all__ = [ + "ScalarField", + "VectorField", + "PatchHierarchy", +] + + +def hierarchy_from( + simulator=None, qty=None, pop="", h5_filename=None, times=None, hier=None, **kwargs +): + 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 times is not None: + times = listify(times) + + 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, times, hier, **kwargs) + + 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..e934bfda5 --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/fromh5.py @@ -0,0 +1,320 @@ +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 .hierarchy import format_timestamp +from ...core.box import Box +from ...core.phare_utilities import ( + refinement_ratio, + none_iterable, + all_iterables, +) +from ...core.gridlayout import GridLayout +from .hierarchy_utils import field_qties +import h5py +from pathlib import Path +from pyphare.core.phare_utilities import listify + + +h5_time_grp_key = "t" +particle_files_patterns = ("domain", "patchGhost", "levelGhost") + + +def get_all_available_quantities_from_h5(filepath, time=0, exclude=["tags"]): + time = format_timestamp(time) + hier = None + path = Path(filepath) + for h5 in path.glob("*.h5"): + if h5.parent == path and h5.stem not in exclude: + hier = hierarchy_fromh5(str(h5), time, hier) + return hier + + +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, as_float=True): + f = h5py.File(filepath, "r") + if as_float: + times = np.array(sorted([float(s) for s in list(f[h5_time_grp_key].keys())])) + else: + times = list(f[h5_time_grp_key].keys()) + f.close() + return times + + +def create_from_all_times(time, hier): + return time is 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): + 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 patch_levels_from_h5(h5f, time, selection_box=None): + """ + creates a dictionary of PatchLevels from a given time in a h5 file + {ilvl: PatchLevel} + """ + + root_cell_width = h5f.attrs["cell_width"] + interp_order = h5f.attrs["interpOrder"] + basename = os.path.basename(h5f.filename) + + patch_levels = {} + + for lvl_key, lvl in h5f[h5_time_grp_key][time].items(): + 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 + + +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 + + h5f = h5py.File(filepath, "r") + selection_box = kwargs.get("selection_box", None) + if hier.has_time(time): + add_data_from_h5(hier, filepath, time) + + 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 + + +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") + + h5f = h5py.File(filepath, "r") + + # force using the hierarchy selection box at that time if existing + 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): + patch.patch_datas.update(patch_levels[ilvl].patches[ip].patch_datas) + + return hier + + +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 + # 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 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[it]) + patch_levels_per_time.append(patch_levels) + + 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_per_time, + domain_box, + refinement_ratio, + times, + h5f, + selection_box=selection_box, + ) + + return hier + + +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_times(time, hier): + time = listify(time) + return new_from_h5(h5_filename, time, **kwargs) + + if create_from_all_times(time, hier): + times = get_times_from_h5(h5_filename) + return new_from_h5(h5_filename, times, **kwargs) + + if load_one_time(time, hier): + time = listify(time) + assert len(time) == 1 + return add_time_from_h5(hier, h5_filename, time[0], **kwargs) + + if load_all_times(time, hier): + 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 + + assert False 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..2e3b62df5 --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -0,0 +1,741 @@ +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 +from ...core.phare_utilities import listify + +import numpy as np +import matplotlib.pyplot as plt + + +def format_timestamp(timestamp): + if isinstance(timestamp, str): + return timestamp + return "{:.10f}".format(timestamp) + + +class PatchHierarchy(object): + """is a collection of patch levels""" + + def __init__( + self, + patch_levels, + domain_box, + refinement_ratio=2, + times=[0.0], + data_files=None, + **kwargs, + ): + if not isinstance(times, (tuple, list)): + times = listify(times) + + if not isinstance(patch_levels, (tuple, list)): + patch_levels = listify(patch_levels) + + 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 = { + format_timestamp(t): box for t, box in zip(times, self.selection_box) + } + assert len(times) == len(self.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( + {format_timestamp(t): pl for t, pl in zip(times, patch_levels)} + ) + + self.domain_box = domain_box + self.refinement_ratio = refinement_ratio + + self._sim = None + + 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} + else: + self.data_files = {} + + self.update() + + def __getitem__(self, qty): + return self.__dict__[qty] + + def update(self): + if len(self.quantities()) > 1: + for qty in self.quantities(): + 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 qty not in self.__dict__: + self.__dict__[qty] = PatchHierarchy( + new_lvls, + self.domain_box, + selection_box=self.domain_box, + times=time, + data_files=self.data_files, + ) + 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: + return self._sim + + # 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( + first_file["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 "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[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 add_time(self, time, patch_level, h5file, selection_box=None): + formated_time = format_timestamp(time) + + self.time_hier[format_timestamp(time)] = patch_level + if selection_box is not None: + self.selection_box[formated_time] = selection_box + + self.data_files[h5file.filename] = h5file + self.update() + + 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(): + print("WARNING - 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 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(): + 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): + 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 has_time(self, time): + return 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()), 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)) + 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..2efadfc3d --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py @@ -0,0 +1,610 @@ +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 + +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): + 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 + + +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]] + return hierarchy + + +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 + """ + + 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_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) + return PatchHierarchy( + patch_levels_per_time, + domain_box, + refinement_ratio, + times=reference_hier.times(), + ) + + +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.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()) + } + return patch_datas + + +def new_patchdatas_from(compute, patchdatas, layout, id, **kwargs): + new_patch_datas = {} + datas = compute(patchdatas, patch_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, t, **kwargs): + reference_hier = hierarchies[0] + new_patches = [] + patch_nbr = len(reference_hier.level(ilvl, time=t).patches) + for ip in range(patch_nbr): + current_patch = reference_hier.level(ilvl, time=t).patches[ip] + layout = current_patch.layout + patch_datas = extract_patchdatas(hierarchies, ilvl, t, 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 + + +def compute_rename(patch_datas, **kwargs): + new_names = kwargs["new_names"] + pd_attrs = [] + + for new_name, pd_name in zip(new_names, patch_datas): + pd_attrs.append( + { + "name": new_name, + "data": patch_datas[pd_name].dataset, + "centering": patch_datas[pd_name].centerings, + } + ) + + return tuple(pd_attrs) + + +def rename(hierarchy, names): + return compute_hier_from(compute_rename, hierarchy, new_names=names) + + +def _compute_mul(patch_datas, **kwargs): + names = kwargs["names"] + + # multiplication of a scalarField or VectorField by a scalar + if "other" in kwargs: + other = kwargs["other"] + pd_attrs = [] + + for name, pd_name in zip(names, patch_datas): + pd_attrs.append( + { + "name": name, + "data": other * patch_datas[pd_name].dataset[:], + "centering": patch_datas[pd_name].centerings, + } + ) + # multiplication of 2 scalarField + elif "self_value" in patch_datas: + dset_value = ( + patch_datas["self_value"].dataset[:] * patch_datas["other_value"].dataset[:] + ) + pd_attrs = ( + { + "name": "value", + "data": dset_value, + "centering": patch_datas["self_value"].centerings, + }, + ) + + return tuple(pd_attrs) + + +def _compute_add(patch_datas, **kwargs): + ref_name = next(iter(patch_datas.keys())) + + dset_x = patch_datas["self_x"].dataset[:] + patch_datas["other_x"].dataset[:] + dset_y = patch_datas["self_y"].dataset[:] + patch_datas["other_y"].dataset[:] + dset_z = patch_datas["self_z"].dataset[:] + patch_datas["other_z"].dataset[:] + + return ( + {"name": "x", "data": dset_x, "centering": patch_datas[ref_name].centerings}, + {"name": "y", "data": dset_y, "centering": patch_datas[ref_name].centerings}, + {"name": "z", "data": dset_z, "centering": patch_datas[ref_name].centerings}, + ) + + +def _compute_sub(patch_datas, **kwargs): + ref_name = next(iter(patch_datas.keys())) + + dset_x = patch_datas["self_x"].dataset[:] - patch_datas["other_x"].dataset[:] + dset_y = patch_datas["self_y"].dataset[:] - patch_datas["other_y"].dataset[:] + dset_z = patch_datas["self_z"].dataset[:] - patch_datas["other_z"].dataset[:] + + return ( + {"name": "x", "data": dset_x, "centering": patch_datas[ref_name].centerings}, + {"name": "y", "data": dset_y, "centering": patch_datas[ref_name].centerings}, + {"name": "z", "data": dset_z, "centering": patch_datas[ref_name].centerings}, + ) + + +def _compute_neg(patch_datas, **kwargs): + names = kwargs["new_names"] + pd_attrs = [] + + for name in names: + pd_attrs.append( + { + "name": name, + "data": -patch_datas[name].dataset[:], + "centering": patch_datas[name].centerings, + } + ) + + return tuple(pd_attrs) + + +def _compute_truediv(patch_datas, **kwargs): + names = kwargs["res_names"] + pd_attrs = [] + + # the denominator is a scalar field which name is "value" + # hence the associated patchdata has to be removed from the list + # of patchdata from the vectorField of the numerator + left_ops = {k: v for k, v in patch_datas.items() if k != "value"} + right_op = patch_datas["value"] + for name, left_op in zip(names, left_ops.values()): + pd_attrs.append( + { + "name": name, + "data": left_op.dataset / right_op.dataset, + "centering": patch_datas[name].centerings, + } + ) + + return tuple(pd_attrs) + + +def _compute_scalardiv(patch_datas, **kwargs): + names = kwargs["res_names"] + scalar = kwargs["scalar"] + pd_attrs = [] + + left_ops = {k: v for k, v in patch_datas.items()} + for name, left_op in zip(names, left_ops.values()): + pd_attrs.append( + { + "name": name, + "data": left_op.dataset / scalar, + "centering": patch_datas[name].centerings, + } + ) + + return tuple(pd_attrs) + + +from dataclasses import dataclass + + +@dataclass +class EqualityReport: + ok: bool + reason: str + + def __bool__(self): + return self.ok + + +def hierarchy_compare(this, that): + if not isinstance(this, PatchHierarchy) or not isinstance(that, PatchHierarchy): + return EqualityReport(False, "class type mismatch") + + if this.ndim != that.ndim or this.domain_box != that.domain_box: + return EqualityReport(False, "dimensional mismatch") + + if this.time_hier.keys() != that.time_hier.keys(): + return EqualityReport(False, "timesteps mismatch") + + for tidx in this.times(): + patch_levels_ref = this.time_hier[tidx] + patch_levels_cmp = that.time_hier[tidx] + + if patch_levels_ref.keys() != patch_levels_cmp.keys(): + return EqualityReport(False, "levels mismatch") + + for level_idx in patch_levels_cmp.keys(): + patch_level_ref = patch_levels_ref[level_idx] + patch_level_cmp = patch_levels_cmp[level_idx] + + for patch_idx in range(len(patch_level_cmp.patches)): + patch_ref = patch_level_ref.patches[patch_idx] + patch_cmp = patch_level_cmp.patches[patch_idx] + + if patch_ref.patch_datas.keys() != patch_cmp.patch_datas.keys(): + print(list(patch_ref.patch_datas.keys())) + print(list(patch_cmp.patch_datas.keys())) + return EqualityReport(False, "data keys mismatch") + + for patch_data_key in patch_ref.patch_datas.keys(): + patch_data_ref = patch_ref.patch_datas[patch_data_key] + patch_data_cmp = patch_cmp.patch_datas[patch_data_key] + + if patch_data_cmp != patch_data_ref: + return EqualityReport( + False, + "data mismatch: " + + type(patch_data_cmp).__name__ + + " " + + type(patch_data_ref).__name__, + ) + + return EqualityReport(True, "OK") diff --git a/pyphare/pyphare/pharesee/hierarchy/patch.py b/pyphare/pyphare/pharesee/hierarchy/patch.py new file mode 100644 index 000000000..0ac18f3d8 --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/patch.py @@ -0,0 +1,69 @@ +# + + +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 __getitem__(self, key): + return self.patch_datas[key] + + 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..717cd10ac --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/patchdata.py @@ -0,0 +1,223 @@ +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 __eq__(self, that): + return self.field_name == that.field_name and self.dataset[:] == that.dataset[:] + + 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_or_slice): + if isinstance(box_or_slice, slice): + return self.dataset[box_or_slice] + return self.select(box_or_slice) + + 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.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() + + def __eq__(self, that): + return self.name == that.name and self.dataset == that.dataset 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/hierarchy/scalarfield.py b/pyphare/pyphare/pharesee/hierarchy/scalarfield.py new file mode 100644 index 000000000..f8691098d --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/scalarfield.py @@ -0,0 +1,212 @@ +from .hierarchy import PatchHierarchy +from .hierarchy_utils import compute_hier_from, compute_rename, rename, _compute_neg + + +class ScalarField(PatchHierarchy): + def __init__(self, hier): + renamed_hier = compute_hier_from(compute_rename, hier, new_names=("value",)) + patch_levels = renamed_hier.patch_levels + domain_box = renamed_hier.domain_box + refinement_ratio = renamed_hier.refinement_ratio + data_files = renamed_hier.data_files + + super().__init__( + patch_levels, domain_box, refinement_ratio, renamed_hier.times(), data_files + ) + + def __add__(self, other): + assert isinstance(other, (ScalarField, int, float)) + h_self = rename(self, ["self_value"]) + + if isinstance(other, ScalarField): + h_other = rename(other, ["other_value"]) + h = compute_hier_from( + self._compute_add, + (h_self, h_other), + ) + elif isinstance(other, (int, float)): + h = compute_hier_from(self._compute_add, (h_self,), other=other) + else: + raise RuntimeError("right operand not supported") + + return ScalarField(h) + + def __radd__(self, other): + return self.__add__(other) + + def __sub__(self, other): + assert isinstance(other, (ScalarField, int, float)) + h_self = rename(self, ["self_value"]) + + if isinstance(other, ScalarField): + h_other = rename(other, ["other_value"]) + h = compute_hier_from( + self._compute_sub, + (h_self, h_other), + ) + elif isinstance(other, (int, float)): + h = compute_hier_from(self._compute_sub, (h_self,), other=other) + else: + raise RuntimeError("right operand not supported") + + return ScalarField(h) + + def __mul__(self, other): + assert isinstance(other, (ScalarField, int, float)) + h_self = rename(self, ["self_value"]) + + if isinstance(other, ScalarField): + h_other = rename(other, ["other_value"]) + h = compute_hier_from(self._compute_mul, (h_self, h_other)) + elif isinstance(other, (int, float)): + h = compute_hier_from(self._compute_mul, (h_self,), other=other) + else: + raise RuntimeError("right operand not supported") + + return ScalarField(h) + + def __rmul__(self, other): + return self.__mul__(other) + + def __truediv__(self, other): + assert isinstance(other, (ScalarField, int, float)) + h_self = rename(self, ["self_value"]) + + if isinstance(other, ScalarField): + h_other = rename(other, ["other_value"]) + h = compute_hier_from(self._compute_truediv, (h_self, h_other)) + elif isinstance(other, (int, float)): + h = compute_hier_from(self._compute_truediv, (h_self,), other=other) + else: + raise RuntimeError("right operand not supported") + + return ScalarField(h) + + def __rtruediv__(self, other): + assert isinstance(other, (int, float)) + h_self = rename(self, ["self_value"]) + + h = compute_hier_from(self._compute_rtruediv, (h_self,), other=other) + + return ScalarField(h) + + def _compute_add(self, patch_datas, **kwargs): + ref_name = next(iter(patch_datas.keys())) + + if "other" in kwargs: + other = kwargs["other"] + dset_value = patch_datas["self_value"].dataset[:] + other + else: + dset_value = ( + patch_datas["self_value"].dataset[:] + + patch_datas["other_value"].dataset[:] + ) + + return ( + { + "name": "value", + "data": dset_value, + "centering": patch_datas[ref_name].centerings, + }, + ) + + def _compute_sub(self, patch_datas, **kwargs): + # subtract a scalar from the dataset of a scalarField + if "other" in kwargs: + other = kwargs["other"] + dset_value = patch_datas["self_value"].dataset[:] - other + # substraction of 2 scalarFields + else: + dset_value = ( + patch_datas["self_value"].dataset[:] + - patch_datas["other_value"].dataset[:] + ) + + return ( + { + "name": "value", + "data": dset_value, + "centering": patch_datas["self_value"].centerings, + }, + ) + + def _compute_mul(self, patch_datas, **kwargs): + # multiplication of a scalarField by a scalar + if "other" in kwargs: + other = kwargs["other"] + pd_attrs = [] + + for pd_name in patch_datas: + pd_attrs.append( + { + "name": "value", + "data": other * patch_datas[pd_name].dataset[:], + "centering": patch_datas[pd_name].centerings, + } + ) + # multiplication of 2 scalarField + else: + dset_value = ( + patch_datas["self_value"].dataset[:] + * patch_datas["other_value"].dataset[:] + ) + pd_attrs = ( + { + "name": "value", + "data": dset_value, + "centering": patch_datas["self_value"].centerings, + }, + ) + + return tuple(pd_attrs) + + def _compute_truediv(self, patch_datas, **kwargs): + # multiplication of a scalarField by a scalar + if "other" in kwargs: + other = kwargs["other"] + pd_attrs = [] + + for pd_name in patch_datas: + pd_attrs.append( + { + "name": "value", + "data": patch_datas[pd_name].dataset[:] / other, + "centering": patch_datas[pd_name].centerings, + } + ) + # multiplication of 2 scalarField + else: + dset_value = ( + patch_datas["self_value"].dataset[:] + / patch_datas["other_value"].dataset[:] + ) + pd_attrs = ( + { + "name": "value", + "data": dset_value, + "centering": patch_datas["self_value"].centerings, + }, + ) + + return tuple(pd_attrs) + + def _compute_rtruediv(self, patch_datas, **kwargs): + # Scalar divided by a scalarField + other = kwargs["other"] + pd_attrs = [] + + for pd_name in patch_datas: + pd_attrs.append( + { + "name": "value", + "data": other / patch_datas[pd_name].dataset[:], + "centering": patch_datas[pd_name].centerings, + } + ) + + return tuple(pd_attrs) + + def __neg__(self): + names_self = self.quantities() + h = compute_hier_from(_compute_neg, self, new_names=names_self) + return ScalarField(h) diff --git a/pyphare/pyphare/pharesee/hierarchy/vectorfield.py b/pyphare/pyphare/pharesee/hierarchy/vectorfield.py new file mode 100644 index 000000000..e9afa59ba --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/vectorfield.py @@ -0,0 +1,100 @@ +from .hierarchy import PatchHierarchy +from .hierarchy_utils import ( + compute_hier_from, + compute_rename, + rename, + _compute_mul, + _compute_add, + _compute_sub, + _compute_truediv, + _compute_scalardiv, +) +from .scalarfield import ScalarField + + +class VectorField(PatchHierarchy): + def __init__(self, hier): + renamed_hier = compute_hier_from( + compute_rename, hier, new_names=("x", "y", "z") + ) + patch_levels = renamed_hier.patch_levels + domain_box = renamed_hier.domain_box + refinement_ratio = renamed_hier.refinement_ratio + data_files = renamed_hier.data_files + + self.names = ["x", "y", "z"] + + super().__init__( + patch_levels, domain_box, refinement_ratio, renamed_hier.times(), data_files + ) + + def __mul__(self, other): + assert isinstance(other, (int, float)) + h = compute_hier_from(_compute_mul, self, names=["x", "y", "z"], other=other) + return VectorField(h) + + def __rmul__(self, other): + return self.__mul__(other) + + def __add__(self, other): + names_self_kept = self.quantities() + names_other_kept = other.quantities() + + if isinstance(other, VectorField): + names_self = ["self_x", "self_y", "self_z"] + names_other = ["other_x", "other_y", "other_z"] + else: + raise RuntimeError("type of hierarchy not yet considered") + + h_self = rename(self, names_self) + h_other = rename(other, names_other) + + h = compute_hier_from( + _compute_add, + (h_self, h_other), + ) + + self = rename(h_self, names_self_kept) # needed ? + other = rename(h_other, names_other_kept) + + return VectorField(h) + + def __sub__(self, other): + names_self_kept = self.quantities() + names_other_kept = other.quantities() + + if isinstance(other, VectorField): + names_self = ["self_x", "self_y", "self_z"] + names_other = ["other_x", "other_y", "other_z"] + else: + raise RuntimeError("type of hierarchy not yet considered") + + h_self = rename(self, names_self) + h_other = rename(other, names_other) + + h = compute_hier_from( + _compute_sub, + (h_self, h_other), + ) + + self = rename(h_self, names_self_kept) + other = rename(h_other, names_other_kept) + + return VectorField(h) + + def __truediv__(self, other): + if not isinstance(other, (ScalarField, int, float)): + raise RuntimeError("type of operand not considered") + + if isinstance(other, ScalarField): + return VectorField( + compute_hier_from( + _compute_truediv, (self, other), res_names=("x", "y", "z") + ) + ) + elif isinstance(other, (int, float)): + return VectorField( + compute_hier_from( + _compute_scalardiv, (self,), res_names=("x", "y", "z"), scalar=other + ) + ) diff --git a/pyphare/pyphare/pharesee/particles.py b/pyphare/pyphare/pharesee/particles.py index f831545ed..a9112048d 100644 --- a/pyphare/pyphare/pharesee/particles.py +++ b/pyphare/pyphare/pharesee/particles.py @@ -85,7 +85,7 @@ def __eq__(self, that): all_assert_sorted(self, that) return True except AssertionError as ex: - print(f"particles.py:Particles::eq failed with:", ex) + print(f"particles.py:Particles::eq failed with: {ex}") print_trace() return False diff --git a/pyphare/pyphare/pharesee/plotting.py b/pyphare/pyphare/pharesee/plotting.py index 249d84624..180db4375 100644 --- a/pyphare/pyphare/pharesee/plotting.py +++ b/pyphare/pyphare/pharesee/plotting.py @@ -245,17 +245,15 @@ 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 r = Run(run_path) time = kwargs.get("time", None) dim = r.GetDl("finest", time).shape[0] interp = kwargs.get("interp", "nearest") - domain = r.GetDomainSize() if qty in ["Bx", "By", "Bz"]: file = os.path.join(run_path, "EM_B.h5") diff --git a/pyphare/pyphare/pharesee/run/__init__.py b/pyphare/pyphare/pharesee/run/__init__.py new file mode 100644 index 000000000..4c8aac32f --- /dev/null +++ b/pyphare/pyphare/pharesee/run/__init__.py @@ -0,0 +1,3 @@ +from .run import Run + +__all__ = ["Run"] diff --git a/pyphare/pyphare/pharesee/run/run.py b/pyphare/pyphare/pharesee/run/run.py new file mode 100644 index 000000000..dd35d8e75 --- /dev/null +++ b/pyphare/pyphare/pharesee/run/run.py @@ -0,0 +1,277 @@ +import os +import numpy as np + +from pyphare.pharesee.hierarchy import hierarchy_from +from pyphare.pharesee.hierarchy import ScalarField, VectorField + +from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from +from pyphare.pharesee.hierarchy.hierarchy_utils import flat_finest_field +from pyphare.core.phare_utilities import listify + +from pyphare.logger import getLogger +from .utils import _current1d +from .utils import ( + _compute_to_primal, + _compute_pop_pressure, + _compute_pressure, + _compute_current, + _compute_divB, + _get_rank, + make_interpolator, +) + + +logger = getLogger(__name__) + +quantities_per_file = { + "EM_B": "B", + "EM_E": "E", + "ions_bulkVelocity": "Vi", + "ions_density": "Ni", + "particle_count": "nppc", +} + + +class Run: + def __init__(self, path): + self.path = path + import glob + + self.available_diags = glob.glob(os.path.join(self.path, "*.h5")) + + 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( + h5_filename=os.path.join(self.path, filename), + times=times, + hier=h, + **kwargs, + ) + + 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 + with the coordinates of the finest grid where the interpolator + can be calculated (that is the return of flat_finest_field) + """ + if merged: + domain = self.GetDomainSize() + dl = self.GetDl() + + # assumes all qties in the hierarchy have the same ghost width + # so take the first patch data of the first patch of the first level.... + nbrGhosts = list(hierarchy.level(0).patches[0].patch_datas.values())[ + 0 + ].ghosts_nbr + merged_qties = {} + for qty in hierarchy.quantities(): + data, coords = flat_finest_field(hierarchy, qty, time=time) + merged_qties[qty] = make_interpolator( + data, coords, interp, domain, dl, qty, nbrGhosts + ) + return merged_qties + else: + return hierarchy + + def GetTags(self, time, merged=False, **kwargs): + hier = self._get_hierarchy(time, "tags.h5") + return self._get(hier, time, merged, "nearest") + + def GetB(self, time, merged=False, interp="nearest", all_primal=True, **kwargs): + hier = self._get_hierarchy(time, "EM_B.h5", **kwargs) + if not all_primal: + return self._get(hier, time, merged, interp) + + h = compute_hier_from(_compute_to_primal, hier, x="Bx", y="By", z="Bz") + return VectorField(h) + + def GetE(self, time, merged=False, interp="nearest", all_primal=True, **kwargs): + hier = self._get_hierarchy(time, "EM_E.h5", **kwargs) + if not all_primal: + return self._get(hier, time, merged, interp) + + h = compute_hier_from(_compute_to_primal, hier, x="Ex", y="Ey", z="Ez") + return VectorField(h) + + def GetMassDensity(self, time, merged=False, interp="nearest", **kwargs): + hier = self._get_hierarchy(time, "ions_mass_density.h5", **kwargs) + return ScalarField(self._get(hier, time, merged, interp)) + + def GetNi(self, time, merged=False, interp="nearest", **kwargs): + hier = self._get_hierarchy(time, "ions_density.h5", **kwargs) + return ScalarField(self._get(hier, time, merged, interp)) + + 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 ScalarField(self._get(hier, time, merged, interp)) + + def GetVi(self, time, merged=False, interp="nearest", **kwargs): + hier = self._get_hierarchy(time, "ions_bulkVelocity.h5", **kwargs) + return VectorField(self._get(hier, time, merged, interp)) + + 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 VectorField(self._get(hier, time, merged, interp)) + + 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, **kwargs), + ) + return self._get(P, time, merged, interp) # should later be a TensorField + + def GetPi(self, time, merged=False, interp="nearest", **kwargs): + M = self._get_hierarchy(time, "ions_momentum_tensor.h5", **kwargs) + massDensity = self.GetMassDensity(time, **kwargs) + Vi = self._get_hierarchy(time, "ions_bulkVelocity.h5", **kwargs) + Pi = compute_hier_from(_compute_pressure, (M, massDensity, Vi)) + return self._get(Pi, time, merged, interp) # should later be a TensorField + + def GetPe(self, time, merged=False, interp="nearest", all_primal=True): + hier = self._get_hierarchy(time, "ions_density.h5") + + Te = hier.sim.electrons.closure.Te + + if not all_primal: + return Te * self._get(hier, time, merged, interp) + + h = compute_hier_from(_compute_to_primal, hier, scalar="rho") + return ScalarField(h) * Te + + def GetJ(self, time, merged=False, interp="nearest", all_primal=True, **kwargs): + B = self.GetB(time, all_primal=False, **kwargs) + J = compute_hier_from(_compute_current, B) + if not all_primal: + return self._get(J, time, merged, interp) + h = compute_hier_from(_compute_to_primal, J, x="Jx", y="Jy", z="Jz") + return VectorField(h) + + def GetDivB(self, time, merged=False, interp="nearest", **kwargs): + B = self.GetB(time, all_primal=False, **kwargs) + db = compute_hier_from(_compute_divB, B) + return ScalarField(self._get(db, time, merged, interp)) + + 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, all_primal=False, **kwargs) + ranks = compute_hier_from(_get_rank, B) + return ScalarField(self._get(ranks, time, merged, interp)) + + 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, **kwargs) + return hier + return self._get_hierarchy(time, filename(pop_name), hier=hier, **kwargs) + + def GetParticleCount(self, time, **kwargs): + c = self._get_hierarchy(time, "particle_count.h5", **kwargs) + return c + + def GetMass(self, pop_name, **kwargs): + list_of_qty = ["density", "flux", "domain", "levelGhost", "patchGhost"] + list_of_mass = [] + + import h5py + + for qty in list_of_qty: + file = os.path.join(self.path, "ions_pop_{}_{}.h5".format(pop_name, qty)) + if os.path.isfile(file): + h5_file = h5py.File(file, "r") + list_of_mass.append(h5_file.attrs["pop_mass"]) + + assert all(m == list_of_mass[0] for m in list_of_mass) + + return list_of_mass[0] + + def GetDomainSize(self, **kwargs): + h5_filename = "EM_B.h5" # _____ TODO : could be another file + + import h5py + + data_file = h5py.File(os.path.join(self.path, h5_filename), "r") + + root_cell_width = np.asarray(data_file.attrs["cell_width"]) + + return (data_file.attrs["domain_box"] + 1) * root_cell_width + + def GetDl(self, level="finest", time=None): + """ + gives the ndarray containing the grid sizes at the given time + for the hierarchy defined in the given run, and for the given level + (default is 'finest', but can also be a int) + + :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" + files = self.available_diags + any_file = files[0] + h5_filename = any_file + import h5py + + data_file = h5py.File(os.path.join(self.path, h5_filename), "r") + + if time is None: + time = float(list(data_file[h5_time_grp_key].keys())[0]) + + hier = self._get_hierarchy(time, h5_filename) + + if level == "finest": + level = hier.finest_level(time) + fac = np.power(hier.refinement_ratio, level) + + root_cell_width = np.asarray(data_file.attrs["cell_width"]) + + return root_cell_width / fac + + def all_times(self): + from glob import glob + import os + import h5py + + path = self.path + files = self.available_diags + ts = {} + for file in files: + basename = os.path.basename(file).split(".")[0] + ff = h5py.File(file) + time_keys = ff["t"].keys() + time = np.zeros(len(time_keys)) + for it, t in enumerate(time_keys): + time[it] = float(t) + ts[quantities_per_file[basename]] = time + return ts + + def times(self, qty): + return self.all_times()[qty] diff --git a/pyphare/pyphare/pharesee/run.py b/pyphare/pyphare/pharesee/run/utils.py similarity index 52% rename from pyphare/pyphare/pharesee/run.py rename to pyphare/pyphare/pharesee/run/utils.py index 842027157..e9cb5b57d 100644 --- a/pyphare/pyphare/pharesee/run.py +++ b/pyphare/pyphare/pharesee/run/utils.py @@ -1,15 +1,6 @@ -import os +from pyphare.core.gridlayout import yee_centering import numpy as np -from .hierarchy import ( - compute_hier_from, - flat_finest_field, - hierarchy_from, -) -from pyphare.logger import getLogger - -logger = getLogger(__name__) - def _current1d(by, bz, xby, xbz): # jx = 0 @@ -134,16 +125,236 @@ def _compute_divB(patchdatas, **kwargs): raise RuntimeError("dimension not implemented") -def _get_rank(patchdatas, **kwargs): +def _ppp_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (primal,primal,primal) to (primal,primal,primal) + centering that is the centering of moments on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, _, _ = _inner_slices(nb_ghosts) + + inner_all = tuple([inner] * ndim) + return inner_all, (inner_all,) + + +def _pdd_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (dual,primal,primal) to (primal,primal,primal) + centering that is the centering of Bx on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, inner_shift_left, inner_shift_right = _inner_slices(nb_ghosts) + + if ndim == 1: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_all,) + elif ndim == 2: + inner_all = tuple([inner] * ndim) + return inner_all, ((inner, inner_shift_left), (inner, inner_shift_right)) + else: + raise RuntimeError("dimension not yet implemented") + + +def _dpd_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (dual,primal,primal) to (primal,primal,primal) + centering that is the centering of By on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, inner_shift_left, inner_shift_right = _inner_slices(nb_ghosts) + + if ndim == 1: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_shift_left, inner_shift_right) + elif ndim == 2: + inner_all = tuple([inner] * ndim) + return inner_all, ((inner_shift_left, inner), (inner_shift_right, inner)) + else: + raise RuntimeError("dimension not yet implemented") + + +def _ddp_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (dual,primal,primal) to (primal,primal,primal) + centering that is the centering of Bz on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, inner_shift_left, inner_shift_right = _inner_slices(nb_ghosts) + + if ndim == 1: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_shift_left, inner_shift_right) + elif ndim == 2: + inner_all = tuple([inner] * ndim) + return inner_all, ( + (inner_shift_left, inner_shift_left), + (inner_shift_left, inner_shift_right), + (inner_shift_right, inner_shift_left), + (inner_shift_right, inner_shift_right), + ) + else: + raise RuntimeError("dimension not yet implemented") + + +def _dpp_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (dual,primal,primal) to (primal,primal,primal) + centering that is the centering of Ex on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, inner_shift_left, inner_shift_right = _inner_slices(nb_ghosts) + + if ndim == 1: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_shift_left, inner_shift_right) + elif ndim == 2: + inner_all = tuple([inner] * ndim) + return inner_all, ((inner_shift_left, inner), (inner_shift_right, inner)) + else: + raise RuntimeError("dimension not yet implemented") + + +def _pdp_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (dual,primal,primal) to (primal,primal,primal) + centering that is the centering of Ey on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, inner_shift_left, inner_shift_right = _inner_slices(nb_ghosts) + + if ndim == 1: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_all,) + elif ndim == 2: + inner_all = tuple([inner] * ndim) + return inner_all, ((inner, inner_shift_left), (inner, inner_shift_right)) + else: + raise RuntimeError("dimension not yet implemented") + + +def _ppd_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (dual,primal,primal) to (primal,primal,primal) + centering that is the centering of Ez on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, _, _ = _inner_slices(nb_ghosts) + + if ndim == 1: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_all,) + elif ndim == 2: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_all,) + else: + raise RuntimeError("dimension not yet implemented") + + +slices_to_primal_ = { + "primal_primal_primal": _ppp_to_ppp_domain_slicing, + "primal_dual_dual": _pdd_to_ppp_domain_slicing, + "dual_primal_dual": _dpd_to_ppp_domain_slicing, + "dual_dual_primal": _ddp_to_ppp_domain_slicing, + "dual_primal_primal": _dpp_to_ppp_domain_slicing, + "primal_dual_primal": _pdp_to_ppp_domain_slicing, + "primal_primal_dual": _ppd_to_ppp_domain_slicing, +} + + +def merge_centerings(pdname): + from pyphare.core.gridlayout import directions + + return "_".join([yee_centering[d][pdname] for d in directions]) + + +def slices_to_primal(pdname, **kwargs): + return slices_to_primal_[merge_centerings(pdname)](**kwargs) + + +def _compute_to_primal(patchdatas, patch_id, **kwargs): + """ + datasets have NaN in their ghosts... might need to be properly filled + with their neighbors already properly projected on primal + """ + + reference_name = next(iter(kwargs.values())) + reference_pd = patchdatas[reference_name] + nb_ghosts = reference_pd.layout.nbrGhosts( + reference_pd.layout.interp_order, "primal" + ) + ndim = reference_pd.box.ndim + + centerings = ["primal"] * ndim + + pd_attrs = [] + for name, pd_name in kwargs.items(): + pd = patchdatas[pd_name] + + ds = pd.dataset + + ds_shape = list(ds.shape) + for i in range(ndim): + if pd.centerings[i] == "dual": + ds_shape[i] += 1 + + # should be something else than nan values when the ghosts cells + # will be filled with correct values coming from the neighbors + ds_all_primal = np.full(ds_shape, np.nan) + ds_ = np.zeros(ds_shape) + + # inner is the slice containing the points that are updated + # in the all_primal dataset + # chunks is a tupls of all the slices coming from the initial dataset + # that are needed to calculate the average for the all_primal dataset + inner, chunks = slices_to_primal(pd_name, nb_ghosts=nb_ghosts, ndim=ndim) + + for chunk in chunks: + ds_[inner] = np.add(ds_[inner], ds[chunk] / len(chunks)) + ds_all_primal[inner] = ds_[inner] + + pd_attrs.append({"name": name, "data": ds_all_primal, "centering": centerings}) + + return tuple(pd_attrs) + + +def _inner_slices(nb_ghosts): + inner = slice(nb_ghosts, -nb_ghosts) + inner_shift_left = slice(nb_ghosts - 1, -nb_ghosts) + inner_shift_right = slice(nb_ghosts, -nb_ghosts + 1) + + return inner, inner_shift_left, inner_shift_right + + +def _get_rank(patchdatas, patch_id, **kwargs): """ make a field dataset cell centered coding the MPI rank rank is obtained from patch global id == "rank#local_patch_id" """ from pyphare.core.box import grow - reference_pd = patchdatas["Bx"] # take Bx as a reference, but could be any other + reference_pd = patchdatas["Bx"] # Bx as a ref, but could be any other ndim = reference_pd.box.ndim - pid = kwargs["id"] layout = reference_pd.layout centering = "dual" @@ -154,7 +365,7 @@ def _get_rank(patchdatas, **kwargs): pass elif ndim == 2: - data = np.zeros(shape) + int(pid.strip("p").split("#")[0]) + data = np.zeros(shape) + int(patch_id.strip("p").split("#")[0]) return ({"name": "rank", "data": data, "centering": [centering] * 2},) else: raise RuntimeError("Not Implemented yet") @@ -277,209 +488,3 @@ def make_interpolator(data, coords, interp, domain, dl, qty, nbrGhosts): raise ValueError("make_interpolator is not yet 3d") return interpolator, finest_coords - - -class Run: - def __init__(self, path, single_hier_for_all_quantities=False): - 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 - - def _get_hierarchy(self, time, filename, hier=None): - t = "{:.10f}".format(time) - - def _get_hier(h): - return hierarchy_from( - h5_filename=os.path.join(self.path, filename), time=t, hier=h - ) - - if self.single_hier_for_all_quantities: - self.hier = _get_hier(self.hier) - return self.hier - return _get_hier(hier) - - def _get(self, hierarchy, time, merged, interp): - """ - if merged=True, will return an interpolator and a tuple of 1d arrays - with the coordinates of the finest grid where the interpolator - can be calculated (that is the return of flat_finest_field) - """ - if merged: - domain = self.GetDomainSize() - dl = self.GetDl() - - # assumes all qties in the hierarchy have the same ghost width - # so take the first patch data of the first patch of the first level.... - nbrGhosts = list(hierarchy.level(0).patches[0].patch_datas.values())[ - 0 - ].ghosts_nbr - merged_qties = {} - for qty in hierarchy.quantities(): - data, coords = flat_finest_field(hierarchy, qty, time=time) - merged_qties[qty] = make_interpolator( - data, coords, interp, domain, dl, qty, nbrGhosts - ) - return merged_qties - else: - return hierarchy - - def GetTags(self, time, merged=False): - hier = self._get_hierarchy(time, "tags.h5") - return self._get(hier, time, merged, "nearest") - - def GetB(self, time, merged=False, interp="nearest"): - hier = self._get_hierarchy(time, "EM_B.h5") - return self._get(hier, time, merged, interp) - - def GetE(self, time, merged=False, interp="nearest"): - hier = self._get_hierarchy(time, "EM_E.h5") - return self._get(hier, time, merged, interp) - - def GetMassDensity(self, time, merged=False, interp="nearest"): - hier = self._get_hierarchy(time, "ions_mass_density.h5") - return self._get(hier, time, merged, interp) - - def GetNi(self, time, merged=False, interp="nearest"): - hier = self._get_hierarchy(time, "ions_density.h5") - 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") - return self._get(hier, time, merged, interp) - - def GetVi(self, time, merged=False, interp="nearest"): - hier = self._get_hierarchy(time, "ions_bulkVelocity.h5") - 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") - 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) - P = compute_hier_from( - _compute_pop_pressure, - (M, V, N), - popname=pop_name, - mass=self.GetMass(pop_name), - ) - 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") - 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) - 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) - db = compute_hier_from(_compute_divB, B) - return self._get(db, time, merged, interp) - - def GetRanks(self, time, merged=False, interp="nearest"): - """ - 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) - ranks = compute_hier_from(_get_rank, B) - return self._get(ranks, time, merged, interp) - - def GetParticles(self, time, pop_name, hier=None): - 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) - return hier - return self._get_hierarchy(time, filename(pop_name), hier=hier) - - def GetMass(self, pop_name): - list_of_qty = ["density", "flux", "domain", "levelGhost", "patchGhost"] - list_of_mass = [] - - import h5py - - for qty in list_of_qty: - file = os.path.join(self.path, "ions_pop_{}_{}.h5".format(pop_name, qty)) - if os.path.isfile(file): - h5_file = h5py.File(file, "r") - list_of_mass.append(h5_file.attrs["pop_mass"]) - - assert all(m == list_of_mass[0] for m in list_of_mass) - - return list_of_mass[0] - - def GetDomainSize(self): - h5_filename = "EM_B.h5" # _____ TODO : could be another file - - import h5py - - data_file = h5py.File(os.path.join(self.path, h5_filename), "r") - - root_cell_width = np.asarray(data_file.attrs["cell_width"]) - - return (data_file.attrs["domain_box"] + 1) * root_cell_width - - def GetDl(self, level="finest", time=None): - """ - gives the ndarray containing the grid sizes at the given time - for the hierarchy defined in the given run, and for the given level - (default is 'finest', but can also be a int) - - :param level: the level at which get the associated grid size - :param time: the time because level depends on it - """ - - h5_time_grp_key = "t" - h5_filename = "EM_B.h5" # _____ TODO : could be another file - - import h5py - - data_file = h5py.File(os.path.join(self.path, h5_filename), "r") - - if time is None: - time = float(list(data_file[h5_time_grp_key].keys())[0]) - - hier = self._get_hierarchy(time, h5_filename) - - if level == "finest": - level = hier.finest_level(time) - fac = np.power(hier.refinement_ratio, level) - - root_cell_width = np.asarray(data_file.attrs["cell_width"]) - - return root_cell_width / fac - - def GetAllAvailableQties(self, time=0, pops=[]): - assert self.single_hier_for_all_quantities # can't work otherwise - - def _try(fn, *args, **kwargs): - try: - fn(*args, **kwargs) - except FileNotFoundError: - # 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) - - for pop in pops: - _try(self.GetFlux, time, pop) - _try(self.GetN, time, pop) - - return self.hier diff --git a/pyphare/pyphare/simulator/simulator.py b/pyphare/pyphare/simulator/simulator.py index 3e6229924..27698ed7a 100644 --- a/pyphare/pyphare/simulator/simulator.py +++ b/pyphare/pyphare/simulator/simulator.py @@ -107,10 +107,10 @@ def setup(self): return self except: import sys + import traceback - print( - 'Exception caught in "Simulator.setup()": {}'.format(sys.exc_info()[0]) - ) + print('Exception caught in "Simulator.setup()": {}'.format(sys.exc_info())) + print(traceback.extract_stack()) raise ValueError("Error in Simulator.setup(), see previous error") def initialize(self): @@ -159,7 +159,7 @@ def advance(self, dt=None): except KeyboardInterrupt as e: self._throw(f"KeyboardInterrupt in simulator.py::advance: \n{e}") - if self._auto_dump() and self.post_advance != None: + if self._auto_dump() and self.post_advance is not None: self.post_advance(self.cpp_sim.currentTime()) return self diff --git a/pyphare/pyphare_tests/test_pharesee/CMakeLists.txt b/pyphare/pyphare_tests/test_pharesee/CMakeLists.txt index 133d31b65..04b7c7aed 100644 --- a/pyphare/pyphare_tests/test_pharesee/CMakeLists.txt +++ b/pyphare/pyphare_tests/test_pharesee/CMakeLists.txt @@ -5,5 +5,7 @@ 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}) +add_python3_test(test-pharesee-hierarchy test_hierarchy.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/__init__.py b/pyphare/pyphare_tests/test_pharesee/__init__.py index 5579bb6a6..51949d4b9 100644 --- a/pyphare/pyphare_tests/test_pharesee/__init__.py +++ b/pyphare/pyphare_tests/test_pharesee/__init__.py @@ -1,16 +1,17 @@ import numpy as np import pyphare.core.box as boxm -from pyphare.core.box import Box, nDBox +from pyphare.core.box import Box from pyphare.core.phare_utilities import listify from pyphare.core.gridlayout import GridLayout, yee_element_is_primal 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 @@ -240,9 +241,6 @@ def build_hierarchy(**kwargs): - nbr_cells - origin - - interp_order - - domain_size - - cell_width - refinement_ratio - refinement_boxes """ @@ -253,9 +251,6 @@ def build_hierarchy(**kwargs): dim = len(origin) if dim > 1: assert len(nbr_cells) == dim - interp_order = kwargs["interp_order"] - domain_size = kwargs["domain_size"] - cell_width = kwargs["cell_width"] refinement_ratio = kwargs["refinement_ratio"] domain_box = boxm.Box([0] * dim, nbr_cells - 1) diff --git a/pyphare/pyphare_tests/test_pharesee/test_geometry.py b/pyphare/pyphare_tests/test_pharesee/test_geometry.py index a2153ce1b..6023f429f 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 @@ -299,7 +292,7 @@ def test_periodic_list(self, refinement_boxes, expected): for ref_box_i, ref_box in enumerate(periodic_boxes): for cmp_box in periodic_boxes[ref_box_i + 1 :]: - self.assertTrue(ref_box * cmp_box == None) + self.assertIs(ref_box * cmp_box, None) self.assertEqual(expected[ilvl], periodic_boxes) diff --git a/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py b/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py index 4dad87b1c..29a8249e3 100644 --- a/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py +++ b/pyphare/pyphare_tests/test_pharesee/test_geometry_2d.py @@ -1,23 +1,14 @@ 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, ) -from pyphare.core.gridlayout import GridLayout, yee_element_is_primal from pyphare_tests.test_pharesee import build_hierarchy @@ -123,7 +114,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])) @@ -171,7 +162,6 @@ def test_large_patchoverlaps(self, expected): level_overlaps = hierarchy_overlaps(hierarchy) ilvl = 0 - lvl = hierarchy.level(ilvl) overlap_boxes = [] self.assertEqual(len(expected), len(level_overlaps[ilvl])) @@ -275,7 +265,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"])) @@ -297,7 +287,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) @@ -386,7 +376,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/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py new file mode 100644 index 000000000..cecfaf56f --- /dev/null +++ b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py @@ -0,0 +1,348 @@ +import unittest +from ddt import ddt +import numpy as np + +from pyphare.simulator.simulator import Simulator +from pyphare.pharesee.run import Run +from pyphare.pharesee.hierarchy import PatchHierarchy +from pyphare.pharesee.hierarchy import ScalarField, VectorField +from pyphare.core.operators import dot, cross, sqrt, modulus, grad +from pyphare.core.box import Box + +diag_outputs = "phare_outputs/" +time_step_nbr = 20 +time_step = 0.005 +final_time = time_step * time_step_nbr +dt = 10 * time_step +nt = int(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().astype(np.float32), times)) + + def test_hierarchy_is_refined(self): + r = Run(diag_outputs) + time = 0.0 + B = r.GetB(time) + self.assertEqual(len(B.levels()), B.levelNbr()) + self.assertEqual(len(B.levels()), 2) + self.assertEqual(len(B.levels(time)), 2) + self.assertEqual(len(B.levels(time)), B.levelNbr(time)) + + def test_can_get_nbytes(self): + r = Run(diag_outputs) + time = 0.0 + B = r.GetB(time) + self.assertGreater(B.nbytes(), 0) + + def test_hierarchy_has_patches(self): + r = Run(diag_outputs) + time = 0.0 + B = r.GetB(time) + self.assertGreater(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.x, PatchHierarchy)) + self.assertTrue(isinstance(B.y, PatchHierarchy)) + self.assertTrue(isinstance(B.z, 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.x.plot(plot_patches=True, ax=ax1) + Bpartial.x.plot(plot_patches=True, ax=ax2) + ax2.add_patch( + Rectangle( + box.lower, + box.shape[0], + box.shape[1], + ec="w", + fc="none", + lw=3, + ) + ) + ax2.set_title(f"{self.id()}") + fig.savefig(f"{self.id()}.png") + + self.assertTrue(isinstance(Bpartial, PatchHierarchy)) + self.assertLess(Bpartial.nbrPatches(), B.nbrPatches()) + + def test_scalarfield_quantities(self): + r = Run(diag_outputs) + time = 0.0 + Ni = r.GetNi(time) + Vi = r.GetVi(time) + self.assertTrue(isinstance(Ni, ScalarField)) + self.assertTrue(isinstance(Vi, VectorField)) + + def test_sum_two_scalarfields(self): + r = Run(diag_outputs) + time = 0.0 + Ni = r.GetNi(time) + Pe = r.GetPe(time) + s = Ni + Pe + self.assertTrue(isinstance(s, ScalarField)) + self.assertEqual(s.quantities(), ["value"]) + + def test_sum_with_scalar(self): + r = Run(diag_outputs) + time = 0.0 + Ni = r.GetNi(time) + s1 = Ni + 0.1 + s2 = 0.1 + Ni + for s in (s1, s2): + self.assertTrue(isinstance(s, ScalarField)) + self.assertEqual(s.quantities(), ["value"]) + + def test_scalarfield_difference(self): + r = Run(diag_outputs) + time = 0.0 + Ni = r.GetNi(time) + Pe = r.GetPe(time) + s1 = Ni - Pe + s2 = -Ni + s3 = Ni - 0.1 + for s in (s1, s2, s3): + self.assertTrue(isinstance(s, ScalarField)) + self.assertEqual(s.quantities(), ["value"]) + + def test_scalarfield_product(self): + r = Run(diag_outputs) + time = 0.0 + Ni = r.GetNi(time) + Pe = r.GetPe(time) + s1 = Ni * Pe + s2 = Ni * 0.1 + for s in (s1, s2): + self.assertTrue(isinstance(s, ScalarField)) + self.assertEqual(s.quantities(), ["value"]) + + def test_scalarfield_division(self): + r = Run(diag_outputs) + time = 0.0 + Ni = r.GetNi(time) + Pe = r.GetPe(time) + s1 = Ni / Pe + s2 = Ni / 0.1 + s3 = Ni / Ni + s4 = 0.1 / Ni + for s in (s1, s2, s3, s4): + self.assertTrue(isinstance(s, ScalarField)) + self.assertEqual(s.quantities(), ["value"]) + + def test_scalarfield_sqrt(self): + r = Run(diag_outputs) + time = 0.0 + Ni = r.GetNi(time) + nisqrt = sqrt(Ni) + self.assertTrue(isinstance(nisqrt, ScalarField)) + + def test_vectorfield_dot_product(self): + r = Run(diag_outputs) + time = 0.0 + B = r.GetB(time) + s1 = dot(B, B) + s2 = modulus(B) + for s in (s1, s2): + self.assertTrue(isinstance(s, ScalarField)) + self.assertEqual(s.quantities(), ["value"]) + + def test_vectorfield_binary_ops(self): + r = Run(diag_outputs) + time = 0.0 + B = r.GetB(time) + V = r.GetVi(time) + s1 = B + V + s2 = B - V + # s3 = -B # TODO + s4 = 2 * B + s5 = B * 2 + s6 = B / 10 + for s in (s1, s2, s4, s5, s6): + self.assertTrue(isinstance(s, VectorField)) + self.assertEqual(s.quantities(), ["x", "y", "z"]) + + def test_scalarfield_to_vectorfield_ops(self): + r = Run(diag_outputs) + time = 0.0 + Ni = r.GetNi(time) + s1 = grad(Ni) + for s in (s1,): + self.assertTrue(isinstance(s, VectorField)) + self.assertEqual(s.quantities(), ["x", "y", "z"]) + + def test_vectorfield_to_vectorfield_ops(self): + r = Run(diag_outputs) + time = 0.0 + B = r.GetB(time) + E = r.GetE(time) + s1 = cross(E, B) + for s in (s1,): + self.assertTrue(isinstance(s, VectorField)) + self.assertEqual(s.quantities(), ["x", "y", "z"]) + + +if __name__ == "__main__": + unittest.main() diff --git a/src/amr/tagging/default_hybrid_tagger_strategy.hpp b/src/amr/tagging/default_hybrid_tagger_strategy.hpp index 6ab30fb07..4ee35e845 100644 --- a/src/amr/tagging/default_hybrid_tagger_strategy.hpp +++ b/src/amr/tagging/default_hybrid_tagger_strategy.hpp @@ -6,6 +6,7 @@ #include "core/data/vecfield/vecfield_component.hpp" #include "core/data/ndarray/ndarray_vector.hpp" #include +#include namespace PHARE::amr 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/core/data/gridlayout/test_laplacian.py b/tests/core/data/gridlayout/test_laplacian.py index 91c689475..8bc8c5754 100644 --- a/tests/core/data/gridlayout/test_laplacian.py +++ b/tests/core/data/gridlayout/test_laplacian.py @@ -210,9 +210,7 @@ def test_laplacian_yee2D(path): Jz = np.tensordot(np.sinh(0.2 * x_primal), np.cosh(0.2 * y_primal), axes=0) Jx_x[psi_d_X : pei_d_X + 1, :] = ( - Jx[ - psi_d_X + 1 : pei_d_X + 2 :, - ] + Jx[psi_d_X + 1 : pei_d_X + 2, :] - 2.0 * Jx[psi_d_X : pei_d_X + 1, :] + Jx[psi_d_X - 1 : pei_d_X, :] ) / (tv.meshSize[0] * tv.meshSize[0]) diff --git a/tests/diagnostic/__init__.py b/tests/diagnostic/__init__.py index a4c4b0cb1..49ab8ef18 100644 --- a/tests/diagnostic/__init__.py +++ b/tests/diagnostic/__init__.py @@ -25,12 +25,13 @@ def dump_all_diags(pops=[], flush_every=100, timestamps=None): write_timestamps=timestamps, compute_timestamps=timestamps, ) - - ph.MetaDiagnostics( - quantity="tags", - write_timestamps=timestamps, - compute_timestamps=timestamps, - ) + # commented out because not working propertly at the moment + # but we keep it there for when it does + # ph.MetaDiagnostics( + # quantity="tags", + # write_timestamps=timestamps, + # compute_timestamps=timestamps, + # ) for quantity in ["density", "bulkVelocity", "pressure_tensor"]: ph.FluidDiagnostics( diff --git a/tests/functional/alfven_wave/alfven_wave1d.py b/tests/functional/alfven_wave/alfven_wave1d.py index aa9ee6846..b0a864629 100644 --- a/tests/functional/alfven_wave/alfven_wave1d.py +++ b/tests/functional/alfven_wave/alfven_wave1d.py @@ -2,9 +2,9 @@ import pyphare.pharein as ph from pyphare.simulator.simulator import Simulator -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 pyphare.pharesee.hierarchy import flat_finest_field +from pyphare.pharesee.hierarchy.hierarchy_utils import flat_finest_field from tests.diagnostic import all_timestamps @@ -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/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): diff --git a/tests/simulator/refinement/test_2d_10_core.py b/tests/simulator/refinement/test_2d_10_core.py index 4b8d50280..07345cb92 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 diff --git a/tests/simulator/refinement/test_2d_2_core.py b/tests/simulator/refinement/test_2d_2_core.py index 4daab318a..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 diff --git a/tests/simulator/test_advance.py b/tests/simulator/test_advance.py index 98701a09b..7fedf997f 100644 --- a/tests/simulator/test_advance.py +++ b/tests/simulator/test_advance.py @@ -17,7 +17,9 @@ ) 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 import format_timestamp +from pyphare.pharesee.hierarchy.hierarchy_utils import merge_particles from pyphare.simulator.simulator import Simulator from tests.diagnostic import all_timestamps @@ -356,7 +358,7 @@ def _test_overlapped_particledatas_have_identical_particles( overlaps = hierarchy_overlaps(datahier, coarsest_time) - for ilvl, lvl in datahier.patch_levels.items(): + for ilvl in datahier.levels(): print("testing level {}".format(ilvl)) for overlap in overlaps[ilvl]: pd1, pd2 = overlap["pdatas"] @@ -454,7 +456,7 @@ def _test_field_coarsening_via_subcycles( uniqTimes = set([0]) for step in range(1, finest_level_step_nbr + 1): - checkTime = datahier.format_timestamp(finestTimeStep * step) + checkTime = format_timestamp(finestTimeStep * step) self.assertIn(checkTime, datahier.times()) uniqTimes.add(checkTime) @@ -470,7 +472,7 @@ def _test_field_coarsening_via_subcycles( ) # skip first coarsest step due to issue 400 for step in range(startStep, syncSteps + 1): - checkTime = datahier.format_timestamp(secondFinestTimeStep * step) + checkTime = format_timestamp(secondFinestTimeStep * step) self.assertIn(checkTime, datahier.times()) nLevels = datahier.levelNbr(checkTime) self.assertGreaterEqual(nLevels, 2) @@ -556,7 +558,7 @@ def base_test_field_level_ghosts_via_subcycles_and_coarser_interpolation( def assert_time_in_hier(*ts): for t in ts: - self.assertIn(L0L1_datahier.format_timestamp(t), L0L1_datahier.times()) + self.assertIn(format_timestamp(t), L0L1_datahier.times()) checks = 0 ndim = global_vars.sim.ndim @@ -710,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( diff --git a/tests/simulator/test_diagnostic_timestamps.py b/tests/simulator/test_diagnostic_timestamps.py index 569da3021..c0c34c4ad 100644 --- a/tests/simulator/test_diagnostic_timestamps.py +++ b/tests/simulator/test_diagnostic_timestamps.py @@ -12,7 +12,9 @@ from ddt import data, ddt from pyphare.core.box import Box1D from pyphare.pharein import ElectromagDiagnostics, ElectronModel -from pyphare.pharesee.hierarchy import h5_filename_from, h5_time_grp_key, hierarchy_from +from pyphare.pharesee.hierarchy import hierarchy_from +from pyphare.pharesee.hierarchy.fromh5 import h5_filename_from +from pyphare.pharesee.hierarchy.hierarchy import format_timestamp from pyphare.simulator.simulator import Simulator @@ -176,9 +178,7 @@ def test_hierarchy_timestamp_cadence(self, refinement_boxes): self.assertEqual(len(time_hier_keys), len(timestamps)) for i, timestamp in enumerate(time_hier_keys): - self.assertEqual( - hier.format_timestamp(timestamps[i]), timestamp - ) + self.assertEqual(format_timestamp(timestamps[i]), timestamp) if __name__ == "__main__": 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..9f4ed8ec4 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.level(1).patches] + ) + + self.assertTrue((nbrPatchGhostPatchDatasOnL1 > 0) == has_patch_ghost) def _test_levelghostparticles_have_correct_split_from_coarser_particle( self, datahier diff --git a/tests/simulator/test_load_balancing.py b/tests/simulator/test_load_balancing.py index bf08256c1..c368c4627 100644 --- a/tests/simulator/test_load_balancing.py +++ b/tests/simulator/test_load_balancing.py @@ -5,7 +5,7 @@ import unittest from ddt import data, ddt, unpack import pyphare.pharein as ph -from pyphare.pharesee.hierarchy import hierarchy_compare +from pyphare.pharesee.hierarchy.hierarchy_utils import hierarchy_compare from pyphare.pharesee.particles import single_patch_per_level_per_pop_from from pyphare.simulator.simulator import Simulator, startMPI @@ -171,7 +171,7 @@ def get_time(path, time=0, pop="protons", datahier=None): from pyphare.pharesee.hierarchy import hierarchy_from return hierarchy_from( - h5_filename=path + f"/ions_pop_{pop}_domain.h5", time=time, hier=datahier + h5_filename=path + f"/ions_pop_{pop}_domain.h5", times=time, hier=datahier ) @@ -221,6 +221,7 @@ def test_raises(self, **lbkwargs): ) # does not get here + @unittest.skip("should change with moments") @data( dict(auto=True), # tolerance checks dict(on_init=True, every=0), # on init only @@ -243,6 +244,7 @@ def test_has_balanced(self, **lbkwargs): tend_sdev = np.std(list(time_info(diag_dir, timestamps[-1]).values())) self.assertLess(tend_sdev, t0_sdev * 0.1) # empirical + @unittest.skip("should change with moments") def test_has_not_balanced_as_defaults(self): if mpi_size == 1: # doesn't make sense return @@ -256,6 +258,7 @@ def test_has_not_balanced_as_defaults(self): tend_sdev = np.std(list(time_info(diag_dir, timestamps[-1]).values())) self.assertGreater(tend_sdev, t0_sdev * 0.1) # empirical + @unittest.skip("should change with moments") def test_compare_is_and_is_not_balanced(self): if mpi_size == 1: # doesn't make sense return diff --git a/tests/simulator/test_restarts.py b/tests/simulator/test_restarts.py index f847ac56c..d0e07edb4 100644 --- a/tests/simulator/test_restarts.py +++ b/tests/simulator/test_restarts.py @@ -19,6 +19,8 @@ from tests.simulator import SimulatorTest from tests.diagnostic import dump_all_diags +from pyphare.pharesee.hierarchy.patchdata import ParticleData +from pyphare.pharesee.hierarchy.fromh5 import get_all_available_quantities_from_h5 def permute(dic, expected_num_levels): @@ -155,26 +157,22 @@ 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) for time in timestamps: checks = 0 - run0 = Run(diag_dir0, single_hier_for_all_quantities=True) - run1 = Run(diag_dir1, single_hier_for_all_quantities=True) + run0 = Run(diag_dir0) - datahier0 = run0.GetAllAvailableQties(time, pops) - datahier1 = run1.GetAllAvailableQties(time, pops) + datahier0 = get_all_available_quantities_from_h5(diag_dir0, time) + datahier1 = get_all_available_quantities_from_h5(diag_dir1, time) self.assertEqual( - datahier0.level(0).patches[0].patch_datas.keys(), - datahier1.level(0).patches[0].patch_datas.keys(), - ) - n_quantities_per_patch = len( - datahier0.level(0).patches[0].patch_datas.keys() + set(datahier0.quantities()), + set(datahier1.quantities()), ) self.assertEqual(len(datahier0.levels()), len(datahier1.levels())) @@ -197,7 +195,8 @@ def count_levels_and_patches(qty): for pd_key, pd0 in patch0.patch_datas.items(): pd1 = patch1.patch_datas[pd_key] self.assertNotEqual(id(pd0), id(pd1)) - if "domain" in pd_key: + self.assertEqual(type(pd0), type(pd1)) + if isinstance(pd1, ParticleData): try: self.assertEqual(pd0.dataset, pd1.dataset) except AssertionError: @@ -208,10 +207,11 @@ def count_levels_and_patches(qty): np.testing.assert_equal(pd0.dataset[:], pd1.dataset[:]) checks += 1 - n_levels, n_patches = count_levels_and_patches(run0.GetB(time)) + n_levels, n_patches = count_levels_and_patches( + run0.GetB(time, all_primal=False) + ) self.assertEqual(n_levels, expected_num_levels) self.assertGreaterEqual(n_patches, n_levels) # at least one patch per level - self.assertEqual(checks, n_quantities_per_patch * n_patches) @data( *permute( diff --git a/tests/simulator/test_run.py b/tests/simulator/test_run.py index c895d142c..8480db130 100644 --- a/tests/simulator/test_run.py +++ b/tests/simulator/test_run.py @@ -20,13 +20,14 @@ startMPI() time_step = 0.005 -final_time = .1 +final_time = 0.1 time_step_nbr = int(final_time / time_step) -timestamps = np.arange(0, final_time+.01, 0.05) +timestamps = np.arange(0, final_time + 0.01, 0.05) diag_dir = "phare_outputs/test_run" plot_dir = Path(f"{diag_dir}_plots") plot_dir.mkdir(parents=True, exist_ok=True) + def config(): L = 0.5 @@ -45,7 +46,7 @@ def config(): diag_options={ "format": "phareh5", "options": {"dir": diag_dir, "mode": "overwrite"}, - } + }, ) def density(x, y): @@ -154,7 +155,9 @@ def vthz(x, y): write_timestamps=timestamps, population_name=pop, ) - ph.FluidDiagnostics(quantity="density", write_timestamps=timestamps, population_name=pop) + ph.FluidDiagnostics( + quantity="density", write_timestamps=timestamps, population_name=pop + ) return sim @@ -173,22 +176,21 @@ def plot(diag_dir): vmax=2e-10, ) run.GetRanks(time).plot( - filename=plot_file_for_qty("Ranks", time), - plot_patches=True, + filename=plot_file_for_qty("Ranks", time), plot_patches=True ) run.GetN(time, pop_name="protons").plot( - filename=plot_file_for_qty("N", time), - plot_patches=True, + filename=plot_file_for_qty("N", time), plot_patches=True ) - for c in ["x","y","z"]: - run.GetB(time).plot( + for c in ["x", "y", "z"]: + run.GetB(time, all_primal=False).plot( filename=plot_file_for_qty(f"b{c}", time), qty=f"B{c}", plot_patches=True, ) run.GetJ(time).plot( + all_primal=False, filename=plot_file_for_qty("jz", time), - qty="Jz", + qty="z", plot_patches=True, vmin=-2, vmax=2, @@ -200,7 +202,10 @@ def assert_file_exists_with_size_at_least(file, size=10000): if not path.exists(): raise FileNotFoundError("file not found: " + file) if path.stat().st_size < size: - raise ValueError("file has unexpected size, possibly corrupt or not written properly: " + file) + raise ValueError( + "file has unexpected size, possibly corrupt or not written properly: " + + file + ) class RunTest(SimulatorTest): @@ -223,10 +228,10 @@ def test_run(self): plot(diag_dir) for time in timestamps: - for q in ["divb","Ranks","N","jz"]: + for q in ["divb", "Ranks", "N", "jz"]: assert_file_exists_with_size_at_least(plot_file_for_qty(q, time)) - for c in ["x","y","z"]: + for c in ["x", "y", "z"]: assert_file_exists_with_size_at_least(plot_file_for_qty(f"b{c}", time)) cpp.mpi_barrier() @@ -234,4 +239,5 @@ def test_run(self): if __name__ == "__main__": import unittest + unittest.main() diff --git a/tools/python3/phloping.py b/tools/python3/phloping.py index fd980dfd8..6a80fa65b 100644 --- a/tools/python3/phloping.py +++ b/tools/python3/phloping.py @@ -38,25 +38,25 @@ def _particles(self): """ Extract particle count per timestep per level from phlop logging """ - 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 + from pyphare.pharesee.hierarchy.fromh5 import get_times_from_h5 + + filepath = self.run.path + "/particle_count.h5" + all_times = get_times_from_h5(filepath) + + particles_per_level_per_time_step = {} + pcount_hier = None + seen_levels = [] + for it, time in enumerate(all_times): + pcount_hier = self.run.GetParticleCount(time) + for ilvl, lvl in pcount_hier.levels(time).items(): + if ilvl not in seen_levels: + seen_levels += [ilvl] + pc = sum([p.attrs["particle_count"] for p in lvl.patches]) + if ilvl not in particles_per_level_per_time_step: + particles_per_level_per_time_step[ilvl] = np.zeros( + len(all_times), dtype=int + ) + 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):