diff --git a/.github/workflows/cmake_macos.yml b/.github/workflows/cmake_macos.yml index 06332b4b9..0c4441e0a 100644 --- a/.github/workflows/cmake_macos.yml +++ b/.github/workflows/cmake_macos.yml @@ -75,10 +75,10 @@ jobs: working-directory: ${{github.workspace}}/build run: | cmake $GITHUB_WORKSPACE -DCMAKE_VERBOSE_MAKEFILE:BOOL=ON \ - -DCMAKE_BUILD_TYPE=RelWithDebInfo \ + -DCMAKE_BUILD_TYPE=Debug \ -DENABLE_SAMRAI_TESTS=OFF -DCMAKE_C_COMPILER_LAUNCHER=ccache \ -DCMAKE_CXX_COMPILER_LAUNCHER=ccache -DlowResourceTests=ON \ - -DCMAKE_CXX_FLAGS="-DPHARE_DIAG_DOUBLES=1 " + -DCMAKE_CXX_FLAGS="-O2 -DPHARE_DIAG_DOUBLES=1" - name: Build working-directory: ${{github.workspace}}/build diff --git a/.github/workflows/cmake_ubuntu.yml b/.github/workflows/cmake_ubuntu.yml index 860e761a8..0b0706d87 100644 --- a/.github/workflows/cmake_ubuntu.yml +++ b/.github/workflows/cmake_ubuntu.yml @@ -84,7 +84,7 @@ jobs: -DCMAKE_C_COMPILER_LAUNCHER=ccache \ -DCMAKE_CXX_COMPILER_LAUNCHER=ccache \ -DlowResourceTests=ON -DdevMode=ON -Dbench=ON \ - -DCMAKE_CXX_FLAGS="-O3 -DPHARE_DIAG_DOUBLES=1 " -Dphare_configurator=ON + -DCMAKE_CXX_FLAGS="-O3 -DPHARE_DIAG_DOUBLES=1" -Dphare_configurator=ON - name: Build working-directory: ${{github.workspace}}/build diff --git a/pyphare/pyphare/core/box.py b/pyphare/pyphare/core/box.py index 9b59b3a1c..3b3d81377 100644 --- a/pyphare/pyphare/core/box.py +++ b/pyphare/pyphare/core/box.py @@ -37,6 +37,10 @@ def nCells(self): """returns the number of cells in the box""" return self.shape.prod() + def size(self): + """returns the number of cells in the box""" + return self.nCells() + def __str__(self): return "Box({},{})".format(self.lower.tolist(), self.upper.tolist()) diff --git a/pyphare/pyphare/core/operators.py b/pyphare/pyphare/core/operators.py index 1bd81fb12..fbcfd2dba 100644 --- a/pyphare/pyphare/core/operators.py +++ b/pyphare/pyphare/core/operators.py @@ -2,60 +2,47 @@ 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())) - +def _compute_dot_product(patch0, hinfo, other, **kwargs): + patch1 = other.level(hinfo.ilvl, hinfo.time)[hinfo.patch_idx] + ref_name = next(iter(patch0.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"][:] + patch0["x"][:] * patch1["x"][:] + + patch0["y"][:] * patch1["y"][:] + + patch0["z"][:] * patch1["z"][:] ) - return ( - {"name": "value", "data": dset, "centering": patch_datas[ref_name].centerings}, - ) + return ({"name": "value", "data": patch0[ref_name].copy_as(dset)},) -def _compute_sqrt(patch_datas, **kwargs): - ref_name = next(iter(patch_datas.keys())) +def _compute_sqrt(patch, **kwargs): + ref_name = next(iter(patch.patch_datas.keys())) - dset = np.sqrt(patch_datas["value"][:]) + dset = np.sqrt(patch["value"][:]) - return ( - {"name": "value", "data": dset, "centering": patch_datas[ref_name].centerings}, - ) + return ({"name": "value", "data": patch[ref_name].copy_as(dset)},) -def _compute_cross_product(patch_datas, **kwargs): - ref_name = next(iter(patch_datas.keys())) +def _compute_cross_product(patch0, hinfo, other, **kwargs): + patch1 = other.level(hinfo.ilvl, hinfo.time)[hinfo.patch_idx] + ref_name = next(iter(patch0.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"][:] - ) + dset_x = patch0["y"][:] * patch1["z"][:] - patch0["z"][:] * patch1["y"][:] + dset_y = patch0["z"][:] * patch1["x"][:] - patch0["x"][:] * patch1["z"][:] + dset_z = patch0["x"][:] * patch1["y"][:] - patch0["y"][:] * patch1["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}, + {"name": "x", "data": patch0[ref_name].copy_as(dset_x)}, + {"name": "y", "data": patch0[ref_name].copy_as(dset_y)}, + {"name": "z", "data": patch0[ref_name].copy_as(dset_z)}, ) -def _compute_grad(patch_data, **kwargs): - ndim = patch_data["value"].box.ndim +def _compute_grad(patch, **kwargs): + ndim = patch["value"].box.ndim nb_ghosts = kwargs["nb_ghosts"] - ds = patch_data["value"].dataset + ds = patch["value"].dataset ds_shape = list(ds.shape) @@ -74,57 +61,28 @@ def _compute_grad(patch_data, **kwargs): 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}, + {"name": "x", "data": patch["value"].copy_as(ds_x)}, + {"name": "y", "data": patch["value"].copy_as(ds_y)}, + {"name": "z", "data": patch["value"].copy_as(ds_z)}, ) -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: +def dot(hier0, hier1, **kwargs): + if not all(isinstance(hier, VectorField) for hier in [hier0, hier1]): raise RuntimeError("type of hierarchy not yet considered") - hl = rename(hier_left, names_left) - hr = rename(hier_right, names_right) + return ScalarField(compute_hier_from(_compute_dot_product, hier0, other=hier1)) - 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: +def cross(hier0, hier1, **kwargs): + if not all(isinstance(hier, VectorField) for hier in [hier0, hier1]): 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) + return VectorField(compute_hier_from(_compute_cross_product, hier0, other=hier1)) def sqrt(hier, **kwargs): - h = compute_hier_from( - _compute_sqrt, - hier, - ) - - return ScalarField(h) + return ScalarField(compute_hier_from(_compute_sqrt, hier)) def modulus(hier): diff --git a/pyphare/pyphare/core/phare_utilities.py b/pyphare/pyphare/core/phare_utilities.py index 29c28bb92..f12fe82d8 100644 --- a/pyphare/pyphare/core/phare_utilities.py +++ b/pyphare/pyphare/core/phare_utilities.py @@ -122,6 +122,18 @@ def __ge__(self, other): return fp_gtr_equal(self.fp, other.fp, self.atol) +class EqualityCheck: + def __init__(self, eq, msg=""): + self.eq = eq + self.msg = msg + + def __bool__(self): + return self.eq + + def __repr__(self): + return self.msg + + def is_fp32(item): if is_nd_array(item): return item.dtype == np.single diff --git a/pyphare/pyphare/pharein/__init__.py b/pyphare/pyphare/pharein/__init__.py index 5a2514275..65955f1aa 100644 --- a/pyphare/pyphare/pharein/__init__.py +++ b/pyphare/pyphare/pharein/__init__.py @@ -328,6 +328,9 @@ def as_paths(rb): ) if len(simulation.diagnostics) > 0: + if simulation.diag_options is not None and "format" in simulation.diag_options: + add_string(diag_path + "format", simulation.diag_options["format"]) + if simulation.diag_options is not None and "options" in simulation.diag_options: add_string( diag_path + "filePath", simulation.diag_options["options"]["dir"] diff --git a/pyphare/pyphare/pharein/diagnostics.py b/pyphare/pyphare/pharein/diagnostics.py index e11801513..bbffad98a 100644 --- a/pyphare/pyphare/pharein/diagnostics.py +++ b/pyphare/pyphare/pharein/diagnostics.py @@ -62,9 +62,9 @@ def validate_timestamps(clazz, key, **kwargs): timestamps = phare_utilities.np_array_ify(kwargs.get(key, [])) if np.any(timestamps < init_time): - raise RuntimeError( - f"Error: timestamp({sim.time_step_nbr}) cannot be less than simulation.init_time({init_time}))" - ) + timestamps = timestamps[timestamps >= init_time] + print(f"Warning: some timestamps below ({init_time}) are filtered)") + if np.any(timestamps > sim.final_time): raise RuntimeError( f"Error: timestamp({sim.time_step_nbr}) cannot be greater than simulation.final_time({sim.final_time}))" diff --git a/pyphare/pyphare/pharein/simulation.py b/pyphare/pyphare/pharein/simulation.py index 21eeebecb..6997d1527 100644 --- a/pyphare/pyphare/pharein/simulation.py +++ b/pyphare/pyphare/pharein/simulation.py @@ -160,31 +160,35 @@ def check_time(**kwargs): and "time_step" not in kwargs ) + if not any([final_and_dt, final_and_nsteps, nsteps_and_dt]): + raise ValueError( + "Error: Specify either 'final_time' and 'time_step' or 'time_step_nbr' and 'time_step'" + + " or 'final_time' and 'time_step_nbr'" + ) + start_time = kwargs.get("restart_options", {}).get("restart_time", 0) + def _final_time(): + if "final_time" in kwargs: + return kwargs["final_time"] + return start_time + kwargs["time_step"] * kwargs["time_step_nbr"] + + final_time = _final_time() + total_time = final_time - start_time + if final_and_dt: - time_step_nbr = int(kwargs["final_time"] / kwargs["time_step"]) - time_step = kwargs["final_time"] / time_step_nbr + time_step_nbr = int(total_time / kwargs["time_step"]) + time_step = total_time / time_step_nbr elif final_and_nsteps: - time_step = kwargs["final_time"] / kwargs["time_step_nbr"] + time_step = total_time / kwargs["time_step_nbr"] time_step_nbr = kwargs["time_step_nbr"] elif nsteps_and_dt: time_step = kwargs["time_step"] time_step_nbr = kwargs["time_step_nbr"] - else: - raise ValueError( - "Error: Specify either 'final_time' and 'time_step' or 'time_step_nbr' and 'time_step'" - + " or 'final_time' and 'time_step_nbr'" - ) - - return ( - time_step_nbr, - time_step, - kwargs.get("final_time", start_time + time_step * time_step_nbr), - ) + return time_step_nbr, time_step, final_time # ------------------------------------------------------------------------------ @@ -481,7 +485,7 @@ def check_directory(directory, key): # diag_options = {"format":"phareh5", "options": {"dir": "phare_ouputs/"}} def check_diag_options(**kwargs): diag_options = kwargs.get("diag_options", None) - formats = ["phareh5"] + formats = ["phareh5", "pharevtkhdf"] if diag_options is not None and "format" in diag_options: if diag_options["format"] not in formats: raise ValueError("Error - diag_options format is invalid") @@ -727,8 +731,10 @@ def wrapper(simulation_object, **kwargs): kwargs["max_nbr_levels"], ) = check_refinement_boxes(ndim, **kwargs) else: - kwargs["max_nbr_levels"] = kwargs.get("max_nbr_levels", None) - assert kwargs["max_nbr_levels"] is not None # this needs setting otherwise + if "max_nbr_levels" not in kwargs: + print("WARNING, 'max_nbr_levels' is not set, defaulting to 1") + kwargs["max_nbr_levels"] = kwargs.get("max_nbr_levels", 1) + kwargs["refinement_boxes"] = None kwargs["tagging_threshold"] = kwargs.get("tagging_threshold", 0.1) diff --git a/pyphare/pyphare/pharesee/hierarchy/__init__.py b/pyphare/pyphare/pharesee/hierarchy/__init__.py index 318b1a24c..1c0103448 100644 --- a/pyphare/pyphare/pharesee/hierarchy/__init__.py +++ b/pyphare/pyphare/pharesee/hierarchy/__init__.py @@ -1,7 +1,13 @@ +# +# +# + from .scalarfield import ScalarField from .vectorfield import VectorField from .hierarchy import PatchHierarchy -from pyphare.core.phare_utilities import listify + + +from pyphare.core import phare_utilities as phut __all__ = [ "ScalarField", @@ -11,11 +17,19 @@ def hierarchy_from( - simulator=None, qty=None, pop="", h5_filename=None, times=None, hier=None, func=None, **kwargs + simulator=None, + qty=None, + pop="", + h5_filename=None, + times=None, + hier=None, + func=None, + **kwargs, ): from .fromh5 import hierarchy_fromh5 from .fromsim import hierarchy_from_sim from .fromfunc import hierarchy_from_func + from .fromvtkhdf5 import hierarchy_fromvtkhdf """ this function reads an HDF5 PHARE file and returns a PatchHierarchy from @@ -29,14 +43,17 @@ def hierarchy_from( """ if times is not None: - times = listify(times) + times = phut.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 h5_filename.endswith(".h5"): + return hierarchy_fromh5(h5_filename, times, hier, **kwargs) + if h5_filename.endswith(".vtkhdf"): + return hierarchy_fromvtkhdf(h5_filename, times, hier, **kwargs) + raise RuntimeError(f"Unknown h5 file type: {h5_filename}") if simulator is not None and qty is not None: return hierarchy_from_sim(simulator, qty, pop=pop) @@ -44,3 +61,19 @@ def hierarchy_from( return hierarchy_from_func(func, hier, **kwargs) raise ValueError("can't make hierarchy") + + +def all_times_from(h5_filename): + if h5_filename.endswith(".h5"): + from .fromh5 import get_times_from_h5 + + return get_times_from_h5(h5_filename) + if h5_filename.endswith(".vtkhdf"): + from .fromvtkhdf5 import get_times_from_h5 + + return get_times_from_h5(h5_filename) + raise RuntimeError(f"Unknown h5 file type: {h5_filename}") + + +def default_time_from(h5_filename): + return all_times_from(h5_filename)[0] diff --git a/pyphare/pyphare/pharesee/hierarchy/for_vtk/__init__.py b/pyphare/pyphare/pharesee/hierarchy/for_vtk/__init__.py new file mode 100644 index 000000000..9679ede0d --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/for_vtk/__init__.py @@ -0,0 +1,2 @@ +# do not make a python file or module called "vtk" +# this might interfer the actual vtk module diff --git a/pyphare/pyphare/pharesee/hierarchy/for_vtk/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/for_vtk/hierarchy.py new file mode 100644 index 000000000..ac171eef2 --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/for_vtk/hierarchy.py @@ -0,0 +1,11 @@ +# +# +# + + +from pyphare.pharesee import hierarchy + + +class PatchHierarchy(hierarchy.PatchHierarchy): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) diff --git a/pyphare/pyphare/pharesee/hierarchy/for_vtk/patch.py b/pyphare/pyphare/pharesee/hierarchy/for_vtk/patch.py new file mode 100644 index 000000000..a1a703ee5 --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/for_vtk/patch.py @@ -0,0 +1,12 @@ +# +# +# + + +from pyphare.pharesee.hierarchy import patch + + +class Patch(patch.Patch): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + self.attrs = {"TO": "DO"} diff --git a/pyphare/pyphare/pharesee/hierarchy/for_vtk/patchdata.py b/pyphare/pyphare/pharesee/hierarchy/for_vtk/patchdata.py new file mode 100644 index 000000000..c32a0c874 --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/for_vtk/patchdata.py @@ -0,0 +1,79 @@ +# +# +# + + +import pyphare.core.box as boxm +from pyphare.pharesee.hierarchy import patchdata + +from pyphare.core import phare_utilities as phut + + +class VtkFieldDatasetAccessor: + def __init__(self, dataset, cmp_idx, offset, box): + self.box = box + self.cmp_idx = cmp_idx + self.offset = offset + self.dataset = dataset + + def __getitem__(self, slice): + # todo finish slice/box lookup + return self.dataset[:, self.cmp_idx][ + self.offset : self.offset + self.box.size() + ].reshape(self.box.shape, order="F") + + @property + def shape(self): + return self.box.shape + + +class VtkFieldData(patchdata.FieldData): + def __init__(self, layout, name, data, cmp_idx, offset, **kwargs): + import h5py + + data_t = type(data) + if not (data_t is h5py.Dataset or data_t is VtkFieldDatasetAccessor): + raise RuntimeError("VtkFieldData only handles vtkhdf datasets") + + box = layout.box + self.field_box = boxm.Box(box.lower, box.upper + 1) + kwargs["ghosts_nbr"] = [0] * layout.box.ndim + kwargs["centering"] = ["primal"] * layout.box.ndim + data = ( + data + if data_t is VtkFieldDatasetAccessor + else VtkFieldDatasetAccessor(data, cmp_idx, offset, self.field_box) + ) + super().__init__(layout, name, data, **kwargs) + + def compare(self, that, atol=1e-16): + """VTK Diagnostics do not have ghosts values!""" + + try: + that_data = ( + that[:] + if all([that.dataset.shape == self.dataset.shape]) + else that[that.box] + ) + phut.assert_fp_any_all_close(self.dataset[:], that_data, atol=atol) + return True + except AssertionError as e: + return phut.EqualityCheck(False, str(e)) + + def __eq__(self, that): + return self.compare(that) + + def __deepcopy__(self, memo): + no_copy_keys = ["dataset"] # do not copy these things + cpy = phut.deep_copy(self, memo, no_copy_keys) + cpy.dataset = self.dataset + return cpy + + def copy_as(self, data=None, **kwargs): + data = self.dataset if data is None else data + if type(data) is VtkFieldDatasetAccessor: + return type(self)( + self.layout, self.field_name, data, data.cmp_idx, data.offset, **kwargs + ) + # make a normal FieldData + return super().copy_as(data, **kwargs) diff --git a/pyphare/pyphare/pharesee/hierarchy/for_vtk/patchlevel.py b/pyphare/pyphare/pharesee/hierarchy/for_vtk/patchlevel.py new file mode 100644 index 000000000..149ef823f --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/for_vtk/patchlevel.py @@ -0,0 +1,11 @@ +# +# +# + + +from pyphare.pharesee.hierarchy import patchlevel + + +class PatchLevel(patchlevel.PatchLevel): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) diff --git a/pyphare/pyphare/pharesee/hierarchy/fromfunc.py b/pyphare/pyphare/pharesee/hierarchy/fromfunc.py index 616852921..37306c718 100644 --- a/pyphare/pyphare/pharesee/hierarchy/fromfunc.py +++ b/pyphare/pyphare/pharesee/hierarchy/fromfunc.py @@ -1,31 +1,35 @@ -from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from +# +# +# import numpy as np +from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from + def ions_mass_density_func1d(x, **kwargs): - masses = kwargs["masses"] # list of float : the ion pop masses - densities = kwargs["densities"] # list of callable : the ion pop density profiles + masses = kwargs["masses"] # list of float : the ion pop masses + densities = kwargs["densities"] # list of callable : the ion pop density profiles assert len(masses) == len(densities) - funcs = np.zeros((x.size, len(masses))) + funcs = np.zeros((x.size, len(masses))) for i, (mass, density) in enumerate(zip(masses, densities)): - funcs[:,i] = mass*density(x) + funcs[:, i] = mass * density(x) return funcs.sum(axis=1) def ions_charge_density_func1d(x, **kwargs): - charges = kwargs["charges"] # list of float : the ion pop charges - densities = kwargs["densities"] # list of callable : the ion pop density profiles + charges = kwargs["charges"] # list of float : the ion pop charges + densities = kwargs["densities"] # list of callable : the ion pop density profiles assert len(charges) == len(densities) - funcs = np.zeros((x.size, len(charges))) + funcs = np.zeros((x.size, len(charges))) for i, (charge, density) in enumerate(zip(charges, densities)): - funcs[:,i] = charge*density(x) + funcs[:, i] = charge * density(x) return funcs.sum(axis=1) @@ -33,43 +37,45 @@ def ions_charge_density_func1d(x, **kwargs): def hierarchy_from_func1d(func, hier, **kwargs): assert hier.ndim == 1 - def compute_(patch_datas, **kwargs): - ref_name = next(iter(patch_datas.keys())) - x_ = patch_datas[ref_name].x + def compute_(patch, **kwargs): + ref_name = next(iter(patch.patch_datas.keys())) + x_ = patch[ref_name].x - return ( - {"name": "value", "data": func(x_, **kwargs), "centering": patch_datas[ref_name].centerings}, + new_patch_data = patch[ref_name].copy_as( + func(x_, **kwargs), ghosts_nbr=patch[ref_name].ghosts_nbr ) + return ({"name": "value", "data": new_patch_data},) + return compute_hier_from(compute_, hier, **kwargs) def ions_mass_density_func2d(x, y, **kwargs): - masses = kwargs["masses"] # list of float : the ion pop masses - densities = kwargs["densities"] # list of callable : the ion pop density profiles + masses = kwargs["masses"] # list of float : the ion pop masses + densities = kwargs["densities"] # list of callable : the ion pop density profiles yv, xv = np.meshgrid(y, x) assert len(masses) == len(densities) - funcs = np.zeros((x.size, y.size, len(masses))) + funcs = np.zeros((x.size, y.size, len(masses))) for i, (mass, density) in enumerate(zip(masses, densities)): - funcs[:,:,i] = mass*density(xv, yv) + funcs[:, :, i] = mass * density(xv, yv) return funcs.sum(axis=2) def ions_charge_density_func2d(x, y, **kwargs): - charges = kwargs["charges"] # list of float : the ion pop charges - densities = kwargs["densities"] # list of callable : the ion pop density profiles + charges = kwargs["charges"] # list of float : the ion pop charges + densities = kwargs["densities"] # list of callable : the ion pop density profiles yv, xv = np.meshgrid(y, x) assert len(charges) == len(densities) - funcs = np.zeros((x.size, y.size, len(charges))) + funcs = np.zeros((x.size, y.size, len(charges))) for i, (charge, density) in enumerate(zip(charges, densities)): - funcs[:,:,i] = charge*density(xv, yv) + funcs[:, :, i] = charge * density(xv, yv) return funcs.sum(axis=2) @@ -77,15 +83,17 @@ def ions_charge_density_func2d(x, y, **kwargs): def hierarchy_from_func2d(func, hier, **kwargs): assert hier.ndim == 2 - def compute_(patch_datas, **kwargs): - ref_name = next(iter(patch_datas.keys())) - x_ = patch_datas[ref_name].x - y_ = patch_datas[ref_name].y + def compute_(patch, **kwargs): + ref_name = next(iter(patch.patch_datas.keys())) + x_ = patch[ref_name].x + y_ = patch[ref_name].y - return ( - {"name": "value", "data": func(x_, y_, **kwargs), "centering": patch_datas[ref_name].centerings}, + new_patch_data = patch[ref_name].copy_as( + func(x_, y_, **kwargs), ghosts_nbr=patch[ref_name].ghosts_nbr ) + return ({"name": "value", "data": new_patch_data},) + return compute_hier_from(compute_, hier, **kwargs) @@ -94,4 +102,3 @@ def hierarchy_from_func(func, hier, **kwargs): return hierarchy_from_func1d(func, hier, **kwargs) if hier.ndim == 2: return hierarchy_from_func2d(func, hier, **kwargs) - diff --git a/pyphare/pyphare/pharesee/hierarchy/fromh5.py b/pyphare/pyphare/pharesee/hierarchy/fromh5.py index 369ac2379..f567e5bcc 100644 --- a/pyphare/pyphare/pharesee/hierarchy/fromh5.py +++ b/pyphare/pyphare/pharesee/hierarchy/fromh5.py @@ -111,7 +111,7 @@ def add_to_patchdata(patch_datas, h5_patch_grp, basename, layout): if is_pop_fluid_file(basename): pdata_name = pop_name(basename) + "_" + pdata_name - if dataset_name in patch_datas: + if pdata_name in patch_datas: raise ValueError("error - {} already in patchdata".format(dataset_name)) patch_datas[pdata_name] = pdata @@ -131,13 +131,11 @@ def h5_filename_from(diagInfo): def get_times_from_h5(filepath, as_float=True): import h5py # see doc/conventions.md section 2.1.1 - 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: + with h5py.File(filepath, "r") as f: times = list(f[h5_time_grp_key].keys()) - f.close() - return times + if as_float: + times = np.array(sorted([float(s) for s in times])) + return times def create_from_all_times(time, hier): diff --git a/pyphare/pyphare/pharesee/hierarchy/fromvtkhdf5.py b/pyphare/pyphare/pharesee/hierarchy/fromvtkhdf5.py new file mode 100644 index 000000000..3c620379b --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/fromvtkhdf5.py @@ -0,0 +1,431 @@ +import numpy as np +from pathlib import Path + +from .for_vtk.patch import Patch +from .for_vtk.patchlevel import PatchLevel +from .for_vtk.patchdata import VtkFieldData as FieldData +from .for_vtk.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 + +from pyphare.core import phare_utilities as phut +from pyphare.core.phare_utilities import listify + +base_path = "VTKHDF" +level_base_path = base_path + "/Level" +step_level_path = base_path + "/Steps/Level" +h5_time_ds_path = "/VTKHDF/Steps/Values" + +_qty_per_filename = { + "EM_E": "EM_E", + "EM_B": "EM_B", + "ions_bulkVelocity": "bulkVelocity", +} + +_vec_fields = { + "EM_E", + "EM_B", + "bulkVelocity", +} + + +def get_path_from(h5file, string): + keys = [v for v in string.split("/") if v] + node = h5file + for key in keys: + node = node[key] + return node + + +class VtkFile: + def __init__(self, path): + if not Path(path).exists(): + raise ValueError(f"ERROR: VTKHDF file does not exist: {path}") + + import h5py # see doc/conventions.md section 2.1.1 + + self.filepath = path + self.file = h5py.File(path, "r") + self._times = None + self._domain_box = None + self._level_spacing = {} + + def level_spacing(self, ilvl=0): + if ilvl in self._level_spacing: + return self._level_spacing[ilvl] + level_group = self._get_path_from(level_base_path + str(ilvl)) + self._level_spacing[ilvl] = level_group.attrs["Spacing"][:] + return self._level_spacing[ilvl] + + def level_boxes(self, ilvl=0, time=0): + dim = self.dimension + time_idx = self.time_idx(time) + num_boxes = self._get_path_from( + step_level_path + str(ilvl) + "/NumberOfAMRBox" + )[time_idx] + box_offset = self._get_path_from(step_level_path + str(ilvl) + "/AMRBoxOffset")[ + time_idx + ] + amr_box_ds = self._get_path_from(level_base_path + str(ilvl) + "/AMRBox") + boxes = [0] * num_boxes + + for bi, boff in enumerate(range(box_offset, box_offset + num_boxes)): + box_vals = amr_box_ds[boff] + lo = np.zeros(dim) + up = np.zeros(dim) + for i, di in enumerate(range(0, dim * 2, 2)): + lo[i] = box_vals[di + 0] + up[i] = box_vals[di + 1] + boxes[bi] = Box(lo, up) + return boxes + + def max_num_levels(self): + n = 0 + for i in range(11): + if f"Level{i}" not in self.file["VTKHDF"]: + break + n += 1 + return n + + def time_idx(self, time): + times = self.times() + ret = np.where(np.isclose(times, float(time), 1e-10)) + return ret[0][0] + + def _amr_box_offset(self, ilvl, time_idx): + return self.file["VTKHDF"]["Steps"][f"Level{ilvl}"]["AMRBoxOffset"][time_idx] + + def _num_boxes(self, ilvl, time_idx): + return self.file["VTKHDF"]["Steps"][f"Level{ilvl}"]["NumberOfAMRBox"][time_idx] + + def _data_offset(self, ilvl, time_idx): + return self.file["VTKHDF"]["Steps"][f"Level{ilvl}"]["PointDataOffset"]["data"][ + time_idx + ] + + def _data_set(self, ilvl): + return self.file["VTKHDF"][f"Level{ilvl}"]["PointData"]["data"] + + @property + def dimension(self): + return self.file.attrs["dimension"] + + @property + def interp_order(self): + return self.file.attrs["interpOrder"] + + @property + def domain_box(self): + if self._domain_box is None: + dbox_attr = self.file.attrs["domain_box"] + self._domain_box = Box([0] * len(dbox_attr), dbox_attr) + return self._domain_box + + def has_time(self, time): + return self.time_idx(time) is not None + + def times(self, reset=False): + if reset: + self._times = None + if self._times is None: + self._times = self._get_path_from(h5_time_ds_path)[:] + return self._times + + def _get_path_from(self, string): + return get_path_from(self.file, string) + + +class VtkPatchLevelInfo: + def __init__(self, vtk_file, ilvl, time_idx): + self.vtk_file = vtk_file + self.ilvl = ilvl + self.time_idx = time_idx + self.num_boxes = vtk_file._num_boxes(ilvl, self.time_idx) + self.data_offset = vtk_file._data_offset(ilvl, self.time_idx) + self.box_offset = vtk_file._amr_box_offset(ilvl, self.time_idx) + self.rolling_data_offset = self.data_offset + + def boxes(self): + dim = self.vtk_file.dimension + amr_box_ds = self.vtk_file.file["VTKHDF"][f"Level{self.ilvl}"]["AMRBox"] + box_vals_list = amr_box_ds[self.box_offset : self.box_offset + self.num_boxes] + boxes = [0] * self.num_boxes + for bi in range(self.num_boxes): + lo = np.zeros(dim) + up = np.zeros(dim) + box_vals = box_vals_list[bi] + for i, di in enumerate(range(0, dim * 2, 2)): + lo[i] = box_vals[di + 0] + up[i] = box_vals[di + 1] + boxes[bi] = Box(lo, up) + return boxes + + def __deepcopy__(self, memo): + no_copy_keys = ["vtk_file"] # do not copy these things + cpy = phut.deep_copy(self, memo, no_copy_keys) + cpy.vtk_file = self.vtk_file + return cpy + + +def get_all_available_quantities_from_h5(filepath, time=0, exclude=["tags"], hier=None): + time = format_timestamp(time) + path = Path(filepath) + for h5 in path.glob("*.vtkhdf"): + if h5.parent == path and h5.stem not in exclude: + hier = hierarchy_fromvtkhdf(str(h5), time, hier) + return hier + + +def make_layout(box, origin, cell_width, interp_order): + return GridLayout(box, 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 is_particle_file(filename): + return False + + +def pop_name(basename): + return Path(".vtkhdf").stem.split("_")[2] + + +def add_to_patchdata(vtk_file, lvl_info, patch_idx, patch_datas, basename, layout): + """ + adds data in the h5_patch_grp in the given PatchData dict + returns True if valid h5 patch found + """ + + X_TIMES = [4, 2, 1][vtk_file.dimension - 1] # data in file is sometimes duplicated + + if is_particle_file(basename): + raise RuntimeError("Particle diagnostics are not supported under vtkhdf") + + def _do_field(cmp_id, qty): + if qty not in field_qties: + raise RuntimeError( + "invalid dataset name : {} is not in {}".format(qty, field_qties) + ) + pdata_name = field_qties[qty] + pdata = FieldData( + layout, + pdata_name, + vtk_file._data_set(lvl_info.ilvl), + cmp_id, + lvl_info.rolling_data_offset, + ) + + if is_pop_fluid_file(basename): + pdata_name = pop_name(basename) + "_" + pdata_name + if pdata_name in patch_datas: + raise ValueError("error - {} already in patchdata".format(qty)) + patch_datas[pdata_name] = pdata + return pdata.field_box.size() # constant as all primal + + data_size = 0 + qty = _qty_per_filename[basename] + if qty in _vec_fields: + for cmp_id, cmp in enumerate(["x", "y", "z"]): + data_size = _do_field(cmp_id, f"{qty}_{cmp}") + else: # assume scalar field + data_size = _do_field(0, qty) + + lvl_info.rolling_data_offset += data_size * X_TIMES + 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 + ".vtkhdf").replace("/", "_")[1:] + + +def get_times_from_h5(filepath, as_float=True): + import h5py # see doc/conventions.md section 2.1.1 + + with h5py.File(filepath, "r") as f: + ds = get_path_from(f, h5_time_ds_path) + times = list(ds[:]) + if as_float: + times = np.array(sorted([float(s) for s in times])) + 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(vtk_file, time, selection_box=None): + """ + creates a dictionary of PatchLevels from a given time in a h5 file + {ilvl: PatchLevel} + """ + + interp_order = vtk_file.interp_order + basename = Path(vtk_file.file.filename).stem + + patch_levels = {} + h5_patch = "" # todo + time_idx = vtk_file.time_idx(time) + for ilvl in range(vtk_file.max_num_levels()): + lvl_cell_width = vtk_file.level_spacing(ilvl) + vtk_patch_level_info = VtkPatchLevelInfo(vtk_file, ilvl, time_idx) + + patches = [] + + for ipatch, patch_box in enumerate(vtk_patch_level_info.boxes()): + origin = patch_box.lower * vtk_file.level_spacing(ilvl) + + intersect = None + if selection_box is not None: + 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 = selection_box * pos_patch_box + + if intersect is not None or selection_box is None: + patch_datas = {} + layout = make_layout(patch_box, origin, lvl_cell_width, interp_order) + + # currently, there is always data for patch + add_to_patchdata( + vtk_file, + vtk_patch_level_info, + ipatch, + patch_datas, + basename, + layout, + ) + patches.append(Patch(patch_datas, h5_patch, layout=layout)) + + if len(patches): + 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 + + vtk_file = VtkFile(filepath) + 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(vtk_file, time, selection_box=selection_box) + hier.add_time(time, patch_levels, vtk_file.file, 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") + + vtk_file = VtkFile(filepath) + + # 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(vtk_file, 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): + 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 = [] + + vtk_file = VtkFile(filepath) + for it, time in enumerate(times): + if isinstance(time, float): + time = f"{time:.10f}" + patch_levels = patch_levels_from_h5( + vtk_file, time, selection_box=selection_box[it] + ) + patch_levels_per_time.append(patch_levels) + + # 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, + vtk_file.domain_box, + refinement_ratio, + times, + vtk_file.file, + selection_box=selection_box, + ) + + return hier + + +def hierarchy_fromvtkhdf(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): + for t in VtkFile(h5_filename).times: + add_time_from_h5(hier, h5_filename, t, **kwargs) + return hier + + assert False diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py index 724eb6409..bceac00db 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy.py @@ -73,9 +73,6 @@ def __deepcopy__(self, memo): no_copy_keys = ["data_files"] # do not copy these things return deep_copy(self, memo, no_copy_keys) - def __getitem__(self, qty): - return self.__dict__[qty] - def update(self): if len(self.quantities()) > 1: for qty in self.quantities(): @@ -432,9 +429,11 @@ def plot1d(self, **kwargs): if qty is None: qty = pdata_names[0] - nbrGhosts = patch.patch_datas[qty].ghosts_nbr - val = patch.patch_datas[qty][patch.box] - x = patch.patch_datas[qty].x[nbrGhosts[0] : -nbrGhosts[0]] + pd = patch[qty] + nbrGhosts = pd.ghosts_nbr + any_ghosts = any(nbrGhosts) + val = pd[patch.box] if any_ghosts else pd[:] + x = pd.x[nbrGhosts[0] : -nbrGhosts[0]] if nbrGhosts[0] > 0 else pd.x label = "L{level}P{patch}".format(level=lvl_nbr, patch=ip) marker = kwargs.get("marker", "") ls = kwargs.get("ls", "--") diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy_compute.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy_compute.py new file mode 100644 index 000000000..8fdc22f23 --- /dev/null +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy_compute.py @@ -0,0 +1,94 @@ +# +# +# + +import operator +from copy import deepcopy + + +def rename(hierarchy, names): + from .hierarchy_utils import compute_hier_from + + return compute_hier_from(compute_rename, hierarchy, new_names=names) + + +def compute_rename(patch, **kwargs): + new_names = kwargs["new_names"] + pd_attrs = [] + + for new_name, pd_name in zip(new_names, patch.patch_datas): + pd_attrs.append({"name": new_name, "data": patch[pd_name]}) + + return tuple(pd_attrs) + + +def compute_mul(patch, **kwargs): + return _compute_copy_op(patch, operator.__mul__, **kwargs) + + +def compute_add(patch, **kwargs): + return _compute_copy_op(patch, operator.__add__, **kwargs) + + +def compute_sub(patch, **kwargs): + return _compute_copy_op(patch, operator.__sub__, **kwargs) + + +def compute_truediv(patch, **kwargs): + return _compute_copy_op(patch, operator.__truediv__, **kwargs) + + +def compute_rtruediv(patch, **kwargs): + return _compute_copy_rop(patch, operator.__truediv__, **kwargs) + + +def _compute_copy_do(patch_data, λ): + new_patch_data = deepcopy(patch_data) + new_patch_data.dataset = λ(patch_data.dataset[:]) + return new_patch_data + + +def drop_ghosts(patch, **kwargs): + pd_attrs = [] + ghosts_nbr = [0] * patch.box.ndim + for name, pd in patch.patch_datas.items(): + data = pd[patch.box] if any(pd.ghosts_nbr) else pd[:] + pd_attrs.append( + { + "name": name, + "data": pd.copy_as(data, ghosts_nbr=ghosts_nbr), + } + ) + return tuple(pd_attrs) + + +class DataAccessor: + def __init__(self, hinfo, other): + self.hinfo = hinfo + self.other = other + + def __getitem__(self, key): + hinfo = self.hinfo + if issubclass(type(self.other), type(hinfo.hier)): + return self.other.level(hinfo.ilvl, hinfo.time)[hinfo.patch_idx][ + key + ].dataset[:] + return self.other + + +def _compute_copy_op(patch, op, hinfo, other, reverse=False): + def _(a, b): + return op(b, a) if reverse else op(a, b) + + data = DataAccessor(hinfo, other) + return tuple( + { + "name": name, + "data": _compute_copy_do(pd, lambda ds: _(ds, data[name])), + } + for name, pd in patch.patch_datas.items() + ) + + +def _compute_copy_rop(patch, op, hinfo, other): + return _compute_copy_op(patch, op, hinfo, other, reverse=True) diff --git a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py index d3856873d..58ca6c81e 100644 --- a/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py +++ b/pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py @@ -1,6 +1,10 @@ -from dataclasses import dataclass, field -from copy import deepcopy +# +# +# + import numpy as np +from copy import deepcopy +from dataclasses import dataclass, field from typing import Any, List, Tuple @@ -134,6 +138,14 @@ def getPatch(hier, point): return patches +@dataclass +class HierarchyAccessor: + hier: PatchHierarchy + time: float + ilvl: int + patch_idx: int + + def compute_hier_from(compute, hierarchies, **kwargs): """return a new hierarchy using the callback 'compute' on all patchdatas of the given hierarchies @@ -142,6 +154,7 @@ def compute_hier_from(compute, hierarchies, **kwargs): 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 = [] @@ -152,6 +165,7 @@ def compute_hier_from(compute, hierarchies, **kwargs): ilvl, new_patches_from(compute, hierarchies, ilvl, t, **kwargs) ) patch_levels_per_time.append(patch_levels) + return PatchHierarchy( patch_levels_per_time, domain_box, @@ -171,12 +185,11 @@ def extract_patchdatas(hierarchies, ilvl, t, ipatch): return patch_datas -def new_patchdatas_from(compute, patchdatas, layout, id, **kwargs): +def new_patchdatas_from(compute, patch, **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 + for data in compute(patch, **kwargs): + assert len(data.keys()) == 2 + new_patch_datas[data["name"]] = data["data"] return new_patch_datas @@ -184,14 +197,16 @@ def new_patches_from(compute, hierarchies, ilvl, t, **kwargs): reference_hier = hierarchies[0] new_patches = [] ref_patches = reference_hier.level(ilvl, time=t).patches - for ip, current_patch in enumerate(ref_patches): - 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 - ) + for ip, ref_patch in enumerate(ref_patches): + patch = deepcopy(ref_patch) + patch.patch_datas = extract_patchdatas(hierarchies, ilvl, t, ip) + hinfo = HierarchyAccessor(reference_hier, t, ilvl, ip) new_patches.append( - Patch(new_patch_datas, current_patch.id, attrs=current_patch.attrs) + Patch( + new_patchdatas_from(compute, patch, hinfo=hinfo, **kwargs), + patch.id, + attrs=patch.attrs, + ) ) return new_patches @@ -446,141 +461,6 @@ def flat_finest_field_2d(hierarchy, qty, time=None): 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) - - @dataclass class EqualityReport: failed: List[Tuple[str, Any, Any]] = field(default_factory=lambda: []) @@ -596,7 +476,7 @@ def __repr__(self): phut.assert_fp_any_all_close(ref[:], cmp[:], atol=1e-16) except AssertionError as e: print(e) - return self.failed[0][0] + return self.failed[0][0] if self.failed else "==" def __call__(self, reason, ref=None, cmp=None): self.failed.append((reason, ref, cmp)) @@ -639,6 +519,9 @@ def hierarchy_compare(this, that, atol=1e-16): patch_ref = patch_level_ref.patches[patch_idx] patch_cmp = patch_level_cmp.patches[patch_idx] + if patch_ref.box != patch_cmp.box: + return eqr("patch box mismatch", patch_ref.box, patch_cmp.box) + if patch_ref.patch_datas.keys() != patch_cmp.patch_datas.keys(): return eqr("data keys mismatch") @@ -646,8 +529,11 @@ def hierarchy_compare(this, that, atol=1e-16): patch_data_ref = patch_ref.patch_datas[patch_data_key] patch_data_cmp = patch_cmp.patch_datas[patch_data_key] - if not patch_data_cmp.compare(patch_data_ref, atol=atol): + ret = patch_data_ref.compare(patch_data_cmp, atol=atol) + if not bool(ret): msg = f"data mismatch: {type(patch_data_ref).__name__} {patch_data_key}" + if type(ret) is not bool: + msg += "\n" + str(ret) eqr(msg, patch_data_cmp, patch_data_ref) if not eqr: diff --git a/pyphare/pyphare/pharesee/hierarchy/patchdata.py b/pyphare/pyphare/pharesee/hierarchy/patchdata.py index e0727b9a2..1cb43738f 100644 --- a/pyphare/pyphare/pharesee/hierarchy/patchdata.py +++ b/pyphare/pyphare/pharesee/hierarchy/patchdata.py @@ -1,8 +1,14 @@ +# +# +# + + import numpy as np -from ...core import phare_utilities as phut + +from ...core import gridlayout from ...core import box as boxm -from ...core.box import Box +from ...core import phare_utilities as phut class PatchData: @@ -24,7 +30,9 @@ def __init__(self, layout, quantity): def __deepcopy__(self, memo): no_copy_keys = ["dataset"] # do not copy these things - return phut.deep_copy(self, memo, no_copy_keys) + cpy = phut.deep_copy(self, memo, no_copy_keys) + cpy.dataset = self.dataset + return cpy class FieldData(PatchData): @@ -35,38 +43,20 @@ class FieldData(PatchData): @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], - ) + self._x = self.yeeCoordsFor(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], - ) + self._y = self.yeeCoordsFor(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], - ) + self._z = self.yeeCoordsFor(2) return self._z def primal_directions(self): @@ -81,9 +71,11 @@ def __repr__(self): return self.__str__() def compare(self, that, atol=1e-16): - return self.field_name == that.field_name and phut.fp_any_all_close( - self.dataset[:], that.dataset[:], atol=atol - ) + try: + phut.assert_fp_any_all_close(self.dataset[:], that.dataset[:], atol=atol) + except AssertionError as e: + return phut.EqualityCheck(False, str(e)) + return self.field_name == that.field_name def __eq__(self, that): return self.compare(that) @@ -96,7 +88,7 @@ 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 + assert isinstance(box, boxm.Box) and box.ndim == self.box.ndim gbox = self.ghost_box.copy() gbox.upper += self.primal_directions() @@ -134,38 +126,15 @@ def __init__(self, layout, field_name, data, **kwargs): self._y = None self._z = None + self.dl = np.asarray(layout.dl) 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.centerings = self._resolve_centering(**kwargs) + self.ghosts_nbr = self._resolve_ghost_nbr(**kwargs) self.ghost_box = boxm.grow(layout.box, self.ghosts_nbr) - self.size = np.copy(self.ghost_box.shape) self.offset = np.zeros(self.ndim) @@ -191,6 +160,53 @@ def grid(): return tuple(g[select] for g in mesh) return mesh + def copy_as(self, data=None, **kwargs): + data = data if data is not None else self.dataset + name = kwargs.get("name", self.field_name) + kwargs["centering"] = kwargs.get("centering", self.centerings) + kwargs["ghosts_nbr"] = kwargs.get("ghosts_nbr", self.ghosts_nbr) + return FieldData(self.layout, name, data, **kwargs) + + def yeeCoordsFor(self, idx): + return self.layout.yeeCoordsFor( + self.field_name, + gridlayout.directions[idx], + withGhosts=any(self.ghosts_nbr) and self.field_name != "tags", + centering=self.centerings[idx], + ) + + def _resolve_ghost_nbr(self, **kwargs): + layout = self.layout + ghosts_nbr = kwargs.get("ghosts_nbr", np.zeros(self.ndim, dtype=int)) + if "ghosts_nbr" not in kwargs: + if self.field_name != "tags": + for i, centering in enumerate(self.centerings): + ghosts_nbr[i] = layout.nbrGhosts(layout.interp_order, centering) + return phut.np_array_ify(ghosts_nbr, layout.box.ndim) + + def _resolve_centering(self, **kwargs): + field_name = self.field_name + if "centering" in kwargs: + if isinstance(kwargs["centering"], list): + assert len(kwargs["centering"]) == self.ndim + return kwargs["centering"] + else: + if self.ndim != 1: + raise ValueError( + "FieldData invalid dimenion for centering argument, expected list for dim > 1" + ) + return phut.listify(kwargs["centering"]) + + if field_name in self.layout.centering["X"]: + directions = ["X", "Y", "Z"][: self.ndim] # drop unused directions + return [ + self.layout.qtyCentering(field_name, direction) + for direction in directions + ] + raise ValueError( + f"centering not specified and cannot be inferred from field name : {field_name}" + ) + class ParticleData(PatchData): """ diff --git a/pyphare/pyphare/pharesee/hierarchy/patchlevel.py b/pyphare/pyphare/pharesee/hierarchy/patchlevel.py index 72fee0d5d..48bafbcc4 100644 --- a/pyphare/pyphare/pharesee/hierarchy/patchlevel.py +++ b/pyphare/pyphare/pharesee/hierarchy/patchlevel.py @@ -1,3 +1,8 @@ +# +# +# + + class PatchLevel: """is a collection of patches""" @@ -13,3 +18,10 @@ def level_range(self): 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 __getitem__(self, idx): + return self.patches[idx] + + @property + def cell_width(self): + return self.patches[0].layout.dl diff --git a/pyphare/pyphare/pharesee/hierarchy/scalarfield.py b/pyphare/pyphare/pharesee/hierarchy/scalarfield.py index f8691098d..54ee8caec 100644 --- a/pyphare/pyphare/pharesee/hierarchy/scalarfield.py +++ b/pyphare/pyphare/pharesee/hierarchy/scalarfield.py @@ -1,10 +1,11 @@ +from . import hierarchy_compute as hc from .hierarchy import PatchHierarchy -from .hierarchy_utils import compute_hier_from, compute_rename, rename, _compute_neg +from .hierarchy_utils import compute_hier_from class ScalarField(PatchHierarchy): def __init__(self, hier): - renamed_hier = compute_hier_from(compute_rename, hier, new_names=("value",)) + renamed_hier = compute_hier_from(hc.compute_rename, hier, new_names=("value",)) patch_levels = renamed_hier.patch_levels domain_box = renamed_hier.domain_box refinement_ratio = renamed_hier.refinement_ratio @@ -15,198 +16,25 @@ def __init__(self, hier): ) 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) + return ScalarField(compute_hier_from(hc.compute_add, self, other=other)) 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) + return ScalarField(compute_hier_from(hc.compute_sub, self, other=other)) 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) + return ScalarField(compute_hier_from(hc.compute_mul, self, other=other)) 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) + return ScalarField(compute_hier_from(hc.compute_truediv, self, other=other)) 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) + return ScalarField(compute_hier_from(hc.compute_rtruediv, self, other=other)) def __neg__(self): - names_self = self.quantities() - h = compute_hier_from(_compute_neg, self, new_names=names_self) - return ScalarField(h) + return self * -1 diff --git a/pyphare/pyphare/pharesee/hierarchy/vectorfield.py b/pyphare/pyphare/pharesee/hierarchy/vectorfield.py index e9afa59ba..60c9098e5 100644 --- a/pyphare/pyphare/pharesee/hierarchy/vectorfield.py +++ b/pyphare/pyphare/pharesee/hierarchy/vectorfield.py @@ -1,21 +1,12 @@ +from . import hierarchy_compute as hc 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 +from .hierarchy_utils import compute_hier_from class VectorField(PatchHierarchy): def __init__(self, hier): renamed_hier = compute_hier_from( - compute_rename, hier, new_names=("x", "y", "z") + hc.compute_rename, hier, new_names=("x", "y", "z") ) patch_levels = renamed_hier.patch_levels domain_box = renamed_hier.domain_box @@ -29,72 +20,20 @@ def __init__(self, hier): ) 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) + if type(other) is VectorField: + raise ValueError( + "VectorField * VectorField is ambiguous, use pyphare.core.operators.dot or .prod" + ) + return VectorField(compute_hier_from(hc.compute_mul, self, other=other)) 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) + return VectorField(compute_hier_from(hc.compute_add, self, other=other)) 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) + return VectorField(compute_hier_from(hc.compute_sub, self, other=other)) 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 - ) - ) + return VectorField(compute_hier_from(hc.compute_truediv, self, other=other)) diff --git a/pyphare/pyphare/pharesee/phare_vtk/__init__.py b/pyphare/pyphare/pharesee/phare_vtk/__init__.py new file mode 100644 index 000000000..9af411dae --- /dev/null +++ b/pyphare/pyphare/pharesee/phare_vtk/__init__.py @@ -0,0 +1,7 @@ +# +# +# + +from .plot import plot + +__all__ = ["plot"] diff --git a/pyphare/pyphare/pharesee/phare_vtk/base.py b/pyphare/pyphare/pharesee/phare_vtk/base.py new file mode 100644 index 000000000..1450f1e3c --- /dev/null +++ b/pyphare/pyphare/pharesee/phare_vtk/base.py @@ -0,0 +1,110 @@ +# +# +# + +import vtk + + +class PhaseOutput: + def __init__(self, **kwargs): + self.kwargs = {**kwargs} + + def __iter__(self): + return self.kwargs.items().__iter__() + + +@staticmethod +def surface_filter(output, **kwargs): + surface = vtk.vtkDataSetSurfaceFilter() + surface.SetInputConnection(output.GetOutputPort()) + return PhaseOutput(surface=surface) + + +@staticmethod +def composite_data_geometry_filter(output, **kwargs): + geom = vtk.vtkCompositeDataGeometryFilter() + geom.SetInputConnection(output.GetOutputPort()) + return PhaseOutput(geom=geom) + + +def poly_data_mapper(output, **kwargs): + array_name = kwargs.get("array_name", "data") + mapper = vtk.vtkPolyDataMapper() + mapper.SetInputConnection(output.GetOutputPort()) + mapper.SelectColorArray(array_name) + mapper.SetScalarModeToUsePointData() + mapper.SetScalarRange( + output.GetOutput().GetPointData().GetArray(array_name).GetRange(0) + ) + return PhaseOutput(mapper=mapper) + + +def all_times_in(reader): + info = reader.GetOutputInformation(0) + + time_key = vtk.vtkStreamingDemandDrivenPipeline.TIME_STEPS() + + if not info.Has(time_key): + raise RuntimeError("VTK Error: cannot get times from file") + + times = info.Get(time_key) + return times + + +def bounds_in(reader): + amr = reader.GetOutput() + bbox = vtk.vtkBoundingBox() + + for level in range(amr.GetNumberOfLevels()): + for idx in range(amr.GetNumberOfDataSets(level)): + ds = amr.GetDataSet(level, idx) + if ds is not None: + b = [0.0] * 6 + ds.GetBounds(b) + bbox.AddBounds(b) + + global_bounds = [0.0] * 6 + bbox.GetBounds(global_bounds) + return global_bounds + + +class VtkFile: + @staticmethod + def _phases(): + return [surface_filter, composite_data_geometry_filter, poly_data_mapper] + + def __init__(self, filename, time=None, array_name="data", phases=_phases()): + if len(phases) == 0: + raise RuntimeError("Error: Zero phases!") + + self.array_name = array_name + self.reader = vtk.vtkHDFReader() + self.reader.SetFileName(filename) + self.reader.Update() + self.reader.UpdateInformation() + + self.times = all_times_in(self.reader) + self.reader.UpdateTimeStep(time if time else self.times[-1]) + + _in = self.reader + for i in range(len(phases)): + ret = phases[i](_in, array_name=array_name) + for key, val in ret: + if hasattr(val, "Update"): + val.Update() + setattr(self, key, val) + _in = val + + self.mapper = _in + self.bounds = bounds_in(self.reader) + self.spacing = self.reader.GetOutput().GetDataSet(0, 0).GetSpacing() + + +class VtkFieldFile(VtkFile): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + +class VtkTensorFieldFile(VtkFile): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) diff --git a/pyphare/pyphare/pharesee/phare_vtk/plot.py b/pyphare/pyphare/pharesee/phare_vtk/plot.py new file mode 100644 index 000000000..b5a81aec8 --- /dev/null +++ b/pyphare/pyphare/pharesee/phare_vtk/plot.py @@ -0,0 +1,38 @@ +# +# +# + + +import vtk +from .base import VtkTensorFieldFile + + +def plot(vtk_file, out_file="vtk.png"): + if isinstance(vtk_file, str): + vtk_file = VtkTensorFieldFile(vtk_file) + vtk_file.geom.GetOutput().GetPointData().SetActiveScalars("data") + + colors = vtk.vtkNamedColors() + + renderer = vtk.vtkRenderer() + renderer.AddActor(vtk.vtkActor(mapper=vtk_file.mapper)) + renderer.SetBackground(colors.GetColor3d("Silver")) + + ren_win = vtk.vtkRenderWindow() + ren_win.AddRenderer(renderer) + ren_win.SetSize(800, 600) + ren_win.OffScreenRenderingOn() # do not flash image on screen + + iren = vtk.vtkRenderWindowInteractor() + iren.SetRenderWindow(ren_win) + + ren_win.Render() + + w2if = vtk.vtkWindowToImageFilter() + w2if.SetInput(ren_win) + w2if.Update() + + writer = vtk.vtkPNGWriter() + writer.SetFileName(out_file) + writer.SetInputConnection(w2if.GetOutputPort()) + writer.Write() diff --git a/pyphare/pyphare/pharesee/plotting.py b/pyphare/pyphare/pharesee/plotting.py index 8f7614919..cef70607f 100644 --- a/pyphare/pyphare/pharesee/plotting.py +++ b/pyphare/pyphare/pharesee/plotting.py @@ -325,7 +325,7 @@ def finest_field_plot(run_path, qty, **kwargs): ax.set_aspect("equal") divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.08) - cb = plt.colorbar(im, ax=ax, cax=cax) + plt.colorbar(im, ax=ax, cax=cax) else: raise ValueError("finest_field_plot not yet ready for 3d") diff --git a/pyphare/pyphare/pharesee/run/run.py b/pyphare/pyphare/pharesee/run/run.py index ce6ca3794..e94f6ea58 100644 --- a/pyphare/pyphare/pharesee/run/run.py +++ b/pyphare/pyphare/pharesee/run/run.py @@ -1,13 +1,20 @@ +# +# +# + import os import glob -import numpy as np +from pathlib import Path +from pyphare.pharesee.hierarchy import all_times_from +from pyphare.pharesee.hierarchy import default_time_from from pyphare.pharesee.hierarchy import hierarchy_from +from pyphare.pharesee.hierarchy import hierarchy_compute as hc from pyphare.pharesee.hierarchy import ScalarField, VectorField +from pyphare.core import phare_utilities as phut 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 ( @@ -23,67 +30,12 @@ logger = getLogger(__name__) -quantities_per_file = { - "EM_B": "B", - "EM_E": "E", - "ions_bulkVelocity": "Vi", - "ions_charge_density": "Ni", - "particle_count": "nppc", -} - class Run: def __init__(self, path, default_time=None): self.path = path self.default_time_ = default_time - 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(time=time) - - # 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 + self.available_diags = self._available_diags() def GetTags(self, time, merged=False, **kwargs): hier = self._get_hierarchy(time, "tags.h5") @@ -92,41 +44,39 @@ def GetTags(self, time, merged=False, **kwargs): def GetB(self, time, merged=False, interp="nearest", all_primal=True, **kwargs): if merged: all_primal = False - hier = self._get_hierarchy(time, "EM_B.h5", **kwargs) + hier = self._get_hier_for(time, "EM_B", **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) + return VectorField(compute_hier_from(_compute_to_primal, hier)) def GetE(self, time, merged=False, interp="nearest", all_primal=True, **kwargs): if merged: all_primal = False - hier = self._get_hierarchy(time, "EM_E.h5", **kwargs) + hier = self._get_hier_for(time, "EM_E", **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) + return VectorField(compute_hier_from(_compute_to_primal, hier)) def GetMassDensity(self, time, merged=False, interp="nearest", **kwargs): - hier = self._get_hierarchy(time, "ions_mass_density.h5", **kwargs) + hier = self._get_hier_for(time, "ions_mass_density", **kwargs) return ScalarField(self._get(hier, time, merged, interp)) def GetNi(self, time, merged=False, interp="nearest", **kwargs): - hier = self._get_hierarchy(time, "ions_charge_density.h5", **kwargs) - return ScalarField(self._get(hier, time, merged, interp)) + hier = self._get_hier_for(time, "ions_charge_density", **kwargs) + return ScalarField(self._get(hier, time, merged, interp, drop_ghosts=True)) def GetN(self, time, pop_name, merged=False, interp="nearest", **kwargs): - hier = self._get_hierarchy(time, f"ions_pop_{pop_name}_density.h5", **kwargs) + hier = self._get_hier_for(time, f"ions_pop_{pop_name}_density", **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)) + hier = self._get_hier_for(time, "ions_bulkVelocity", **kwargs) + return VectorField(self._get(hier, time, merged, interp, drop_ghosts=True)) def GetFlux(self, time, pop_name, merged=False, interp="nearest", **kwargs): - hier = self._get_hierarchy(time, f"ions_pop_{pop_name}_flux.h5", **kwargs) + hier = self._get_hier_for(time, f"ions_pop_{pop_name}_flux", **kwargs) return VectorField(self._get(hier, time, merged, interp)) def GetPressure(self, time, pop_name, merged=False, interp="nearest", **kwargs): @@ -146,19 +96,16 @@ def GetPressure(self, time, pop_name, merged=False, interp="nearest", **kwargs): 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) + Vi = self._get_hier_for(time, "ions_bulkVelocity", **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_charge_density.h5") - + hier = self._get_hier_for(time, "ions_charge_density") 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") + h = compute_hier_from(hc.drop_ghosts, hier) return ScalarField(h) * Te def GetJ(self, time, merged=False, interp="nearest", all_primal=True, **kwargs): @@ -168,8 +115,8 @@ def GetJ(self, time, merged=False, interp="nearest", all_primal=True, **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) + J = hc.rename(compute_hier_from(_compute_to_primal, J), ["Jx", "Jy", "Jz"]) + return VectorField(J) def GetDivB(self, time, merged=False, interp="nearest", **kwargs): B = self.GetB(time, all_primal=False, **kwargs) @@ -218,12 +165,10 @@ def GetMass(self, pop_name, **kwargs): return list_of_mass[0] def GetDomainSize(self, **kwargs): - import h5py - - data_file = h5py.File(self.available_diags[0], "r") # That is the first file in th available diags - root_cell_width = np.asarray(data_file.attrs["cell_width"]) - - return (data_file.attrs["domain_box"] + 1) * root_cell_width + hier = self._get_any_hierarchy(self.default_time) + root_cell_width = hier.level(0).cell_width + domain_box = hier.domain_box + return (domain_box.upper + 1) * root_cell_width def GetDl(self, level="finest", time=None): """ @@ -235,55 +180,93 @@ def GetDl(self, level="finest", time=None): :param time: the time because level depends on it """ - import h5py - - def _get_time(): - if time: - return time - if self.default_time_: - return self.default_time_ - self.default_time_ = float(list(data_file[h5_time_grp_key].keys())[0]) - return self.default_time_ + time = time if time is not None else self.default_time + hier = self._get_any_hierarchy(time) + level = hier.finest_level(time) if level == "finest" else level + return hier.level(level).cell_width - h5_time_grp_key = "t" - files = self.available_diags - - for h5_filename in files: - data_file = h5py.File(h5_filename, "r") - - time = _get_time() + def all_times(self): + return {Path(file).stem: all_times_from(file) for file in self.available_diags} - try: - hier = self._get_hierarchy(time, h5_filename.split("/")[-1]) + def times(self, qty): + return self.all_times()[qty] - if level == "finest": - level = hier.finest_level(time) - fac = np.power(hier.refinement_ratio, level) + def _available_diags(self): + files = glob.glob(os.path.join(self.path, "*.h5")) + if files: + return files + files = glob.glob(os.path.join(self.path, "*.vtkhdf")) + if files: + return files + raise RuntimeError(f"No HDF5 files found at: {self.path}") + + def _get_hier_for(self, time, qty, **kwargs): + path = Path(self.path) / f"{qty}.h5" + if path.exists(): + return self._get_hierarchy(time, f"{qty}.h5", **kwargs) + path = Path(self.path) / f"{qty}.vtkhdf" + if path.exists(): + return self._get_hierarchy(time, f"{qty}.vtkhdf", **kwargs) + raise RuntimeError(f"No HDF5 file found for: {qty}") + + def _get_any_hierarchy(self, time): + ref_file = Path(self.available_diags[0]).stem + return self._get_hier_for(time, ref_file) - root_cell_width = np.asarray(data_file.attrs["cell_width"]) + def _get_hierarchy(self, times, filename, hier=None, **kwargs): + from pyphare.core.box import Box - return root_cell_width / fac + times = phut.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) - except KeyError: - ... # time may not be avilaable for given quantity + def _get_hier(h): + return hierarchy_from( + h5_filename=os.path.join(self.path, filename), + times=times, + hier=h, + **kwargs, + ) - raise RuntimeError("Unable toGetDl") + return _get_hier(hier) - def all_times(self): - import h5py + # TODO maybe transform that so multiple times can be accepted + def _get(self, hierarchy, time, merged, interp, drop_ghosts=False): + """ + 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(time=time) - 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 - ff.close() - return ts + # 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 ( + compute_hier_from(hc.drop_ghosts, hierarchy) + if drop_ghosts + else hierarchy + ) - def times(self, qty): - return self.all_times()[qty] + @property + def default_time(self): + if self.default_time_ is None: + ref_file = self.available_diags[0] + self.default_time_ = default_time_from(ref_file) + return self.default_time_ diff --git a/pyphare/pyphare/pharesee/run/utils.py b/pyphare/pyphare/pharesee/run/utils.py index d6ffac24c..d37848855 100644 --- a/pyphare/pyphare/pharesee/run/utils.py +++ b/pyphare/pyphare/pharesee/run/utils.py @@ -1,7 +1,13 @@ -from pyphare.core.gridlayout import yee_centering +# +# +# + import numpy as np +from pyphare.core.gridlayout import yee_centering + + def _current1d(by, bz, xby, xbz): # jx = 0 # jy = -dxBz @@ -59,45 +65,45 @@ def _current2d(bx, by, bz, dx, dy): return jx, jy, jz -def _compute_current(patchdatas, **kwargs): - reference_pd = patchdatas["Bx"] # take Bx as a reference, but could be any other +def _compute_current(patch, **kwargs): + reference_pd = patch["Bx"] # take Bx as a reference, but could be any other ndim = reference_pd.box.ndim if ndim == 1: - By = patchdatas["By"].dataset[:] - xby = patchdatas["By"].x - Bz = patchdatas["Bz"].dataset[:] - xbz = patchdatas["Bz"].x - Jy, Jz = _current1d(By, Bz, xby, xbz) + By = patch["By"] + xby = By.x + Bz = patch["Bz"] + xbz = Bz.x + Jy, Jz = _current1d(By[:], Bz[:], xby, xbz) return ( - {"name": "Jy", "data": Jy, "centering": "primal"}, - {"name": "Jz", "data": Jz, "centering": "primal"}, + {"name": "Jy", "data": By.copy_as(Jy, centering="primal")}, + {"name": "Jz", "data": Bz.copy_as(Jz, centering="primal")}, ) elif ndim == 2: - Bx = patchdatas["Bx"].dataset[:] - By = patchdatas["By"].dataset[:] - Bz = patchdatas["Bz"].dataset[:] + Bx = patch["Bx"] + By = patch["By"] + Bz = patch["Bz"] dx, dy = reference_pd.dl + Jx, Jy, Jz = _current2d(Bx[:], By[:], Bz[:], dx, dy) - Jx, Jy, Jz = _current2d(Bx, By, Bz, dx, dy) - - components = ("Jx", "Jy", "Jz") centering = { component: [ reference_pd.layout.centering[direction][component] for direction in ("X", "Y") ] - for component in components + for component in ("Jx", "Jy", "Jz") } return ( - {"name": "Jx", "data": Jx, "centering": centering["Jx"]}, - {"name": "Jy", "data": Jy, "centering": centering["Jy"]}, - {"name": "Jz", "data": Jz, "centering": centering["Jz"]}, + {"name": "Jx", "data": Bx.copy_as(Jx, centering=centering["Jx"])}, + {"name": "Jy", "data": By.copy_as(Jy, centering=centering["Jy"])}, + {"name": "Jz", "data": Bz.copy_as(Jz, centering=centering["Jz"])}, ) + raise RuntimeError("dimension not implemented") + def _divB2D(Bx, By, xBx, yBy): dxbx = (Bx[1:, :] - Bx[:-1, :]) / (xBx[1] - xBx[0]) @@ -105,24 +111,24 @@ def _divB2D(Bx, By, xBx, yBy): return dxbx + dyby -def _compute_divB(patchdatas, **kwargs): - reference_pd = patchdatas["Bx"] # take Bx as a reference, but could be any other +def _compute_divB(patch, **kwargs): + reference_pd = patch["Bx"] # take Bx as a reference, but could be any other ndim = reference_pd.box.ndim if ndim == 1: raise ValueError("divB is 0 by construction in 1D") - elif ndim == 2: - By = patchdatas["By"].dataset[:] - Bx = patchdatas["Bx"].dataset[:] - xBx = patchdatas["Bx"].x - yBy = patchdatas["By"].y - divB = _divB2D(Bx, By, xBx, yBy) + centering = ["dual"] * ndim - return ({"name": "divB", "data": divB, "centering": ["dual", "dual"]},) + if ndim == 2: + By = patch["By"].dataset[:] + Bx = patch["Bx"].dataset[:] + xBx = patch["Bx"].x + yBy = patch["By"].y + divB = reference_pd.copy_as(_divB2D(Bx, By, xBx, yBy), centering=centering) + return ({"name": "divB", "data": divB},) - else: - raise RuntimeError("dimension not implemented") + raise RuntimeError("dimension not implemented") def _ppp_to_ppp_domain_slicing(**kwargs): @@ -154,11 +160,12 @@ def _pdd_to_ppp_domain_slicing(**kwargs): if ndim == 1: inner_all = tuple([inner] * ndim) return inner_all, (inner_all,) - elif ndim == 2: + + if 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") + + raise RuntimeError("dimension not yet implemented") def _dpd_to_ppp_domain_slicing(**kwargs): @@ -292,48 +299,49 @@ def slices_to_primal(pdname, **kwargs): return slices_to_primal_[merge_centerings(pdname)](**kwargs) -def _compute_to_primal(patchdatas, patch_id, **kwargs): +def _compute_to_primal(patch, **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 + ndim = patch.box.ndim pd_attrs = [] - for name, pd_name in kwargs.items(): - pd = patchdatas[pd_name] + for name, ref_pd in patch.patch_datas.items(): + nb_ghosts = ref_pd.layout.nbrGhosts(ref_pd.layout.interp_order, "primal") + ref_ds = ref_pd.dataset - ds = pd.dataset + should_skip = all( # vtkhdf is all primal with no ghosts + [ref_pd.centerings == ["primal"] * ndim, not any(ref_pd.ghosts_nbr)] + ) + if should_skip: + pd_attrs.append({"name": name, "data": ref_pd}) + continue - ds_shape = list(ds.shape) + ds_shape = list(ref_ds.shape) for i in range(ndim): - if pd.centerings[i] == "dual": + if ref_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) + inner, chunks = slices_to_primal(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] + ds_[inner] = np.add(ds_[inner], ref_ds[chunk] / len(chunks)) + + copy_pd = ref_pd.copy_as( + ds_[inner], centering=["primal"] * ndim, ghosts_nbr=[0] * ndim + ) - pd_attrs.append({"name": name, "data": ds_all_primal, "centering": centerings}) + pd_attrs.append({"name": name, "data": copy_pd}) return tuple(pd_attrs) @@ -346,18 +354,18 @@ def _inner_slices(nb_ghosts): return inner, inner_shift_left, inner_shift_right -def _get_rank(patchdatas, patch_id, **kwargs): +def _get_rank(patch, **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"] # Bx as a ref, but could be any other + reference_pd = patch["Bx"] # Bx as a ref, but could be any other ndim = reference_pd.box.ndim layout = reference_pd.layout - centering = "dual" + centering = ["dual"] * ndim nbrGhosts = layout.nbrGhosts(layout.interp_order, centering) shape = grow(reference_pd.box, [nbrGhosts] * 2).shape @@ -365,42 +373,43 @@ def _get_rank(patchdatas, patch_id, **kwargs): pass elif ndim == 2: - data = np.zeros(shape) + int(patch_id.strip("p").split("#")[0]) - return ({"name": "rank", "data": data, "centering": [centering] * 2},) + data = np.zeros(shape) + int(patch.id.strip("p").split("#")[0]) + pd = reference_pd.copy_as(data, centering=centering) + return ({"name": "rank", "data": pd},) else: raise RuntimeError("Not Implemented yet") -def _compute_pressure(patch_datas, **kwargs): - Mxx = patch_datas["Mxx"].dataset[:] - Mxy = patch_datas["Mxy"].dataset[:] - Mxz = patch_datas["Mxz"].dataset[:] - Myy = patch_datas["Myy"].dataset[:] - Myz = patch_datas["Myz"].dataset[:] - Mzz = patch_datas["Mzz"].dataset[:] - massDensity = patch_datas["value"].dataset[:] - Vix = patch_datas["Vx"].dataset[:] - Viy = patch_datas["Vy"].dataset[:] - Viz = patch_datas["Vz"].dataset[:] - - Pxx = Mxx - Vix * Vix * massDensity - Pxy = Mxy - Vix * Viy * massDensity - Pxz = Mxz - Vix * Viz * massDensity - Pyy = Myy - Viy * Viy * massDensity - Pyz = Myz - Viy * Viz * massDensity - Pzz = Mzz - Viz * Viz * massDensity +def _compute_pressure(patch, **kwargs): + Mxx = patch["Mxx"] + Mxy = patch["Mxy"] + Mxz = patch["Mxz"] + Myy = patch["Myy"] + Myz = patch["Myz"] + Mzz = patch["Mzz"] + massDensity = patch["value"][:] + Vix = patch["Vx"][:] + Viy = patch["Vy"][:] + Viz = patch["Vz"][:] + + Pxx = Mxx[:] - Vix * Vix * massDensity + Pxy = Mxy[:] - Vix * Viy * massDensity + Pxz = Mxz[:] - Vix * Viz * massDensity + Pyy = Myy[:] - Viy * Viy * massDensity + Pyz = Myz[:] - Viy * Viz * massDensity + Pzz = Mzz[:] - Viz * Viz * massDensity return ( - {"name": "Pxx", "data": Pxx, "centering": ["primal", "primal"]}, - {"name": "Pxy", "data": Pxy, "centering": ["primal", "primal"]}, - {"name": "Pxz", "data": Pxz, "centering": ["primal", "primal"]}, - {"name": "Pyy", "data": Pyy, "centering": ["primal", "primal"]}, - {"name": "Pyz", "data": Pyz, "centering": ["primal", "primal"]}, - {"name": "Pzz", "data": Pzz, "centering": ["primal", "primal"]}, + {"name": "Pxx", "data": Mxx.copy_as(Pxx)}, + {"name": "Pxy", "data": Mxy.copy_as(Pxy)}, + {"name": "Pxz", "data": Mxz.copy_as(Pxz)}, + {"name": "Pyy", "data": Myy.copy_as(Pyy)}, + {"name": "Pyz", "data": Myz.copy_as(Pyz)}, + {"name": "Pzz", "data": Mzz.copy_as(Pzz)}, ) -def _compute_pop_pressure(patch_datas, **kwargs): +def _compute_pop_pressure(patch, **kwargs): """ computes the pressure tensor for a given population this method is different from _compute_pressure in that: @@ -411,33 +420,33 @@ def _compute_pop_pressure(patch_datas, **kwargs): P = M - F*F/N * mass """ popname = kwargs["popname"] - Mxx = patch_datas[popname + "_Mxx"].dataset[:] - Mxy = patch_datas[popname + "_Mxy"].dataset[:] - Mxz = patch_datas[popname + "_Mxz"].dataset[:] - Myy = patch_datas[popname + "_Myy"].dataset[:] - Myz = patch_datas[popname + "_Myz"].dataset[:] - Mzz = patch_datas[popname + "_Mzz"].dataset[:] - Fx = patch_datas["x"].dataset[:] - Fy = patch_datas["y"].dataset[:] - Fz = patch_datas["z"].dataset[:] - N = patch_datas["value"].dataset[:] + Mxx = patch[popname + "_Mxx"] + Mxy = patch[popname + "_Mxy"] + Mxz = patch[popname + "_Mxz"] + Myy = patch[popname + "_Myy"] + Myz = patch[popname + "_Myz"] + Mzz = patch[popname + "_Mzz"] + Fx = patch["x"][:] + Fy = patch["y"][:] + Fz = patch["z"][:] + N = patch["value"][:] mass = kwargs["mass"] - Pxx = Mxx - Fx * Fx * mass / N - Pxy = Mxy - Fx * Fy * mass / N - Pxz = Mxz - Fx * Fz * mass / N - Pyy = Myy - Fy * Fy * mass / N - Pyz = Myz - Fy * Fz * mass / N - Pzz = Mzz - Fz * Fz * mass / N + Pxx = Mxx[:] - Fx * Fx * mass / N + Pxy = Mxy[:] - Fx * Fy * mass / N + Pxz = Mxz[:] - Fx * Fz * mass / N + Pyy = Myy[:] - Fy * Fy * mass / N + Pyz = Myz[:] - Fy * Fz * mass / N + Pzz = Mzz[:] - Fz * Fz * mass / N return ( - {"name": popname + "_Pxx", "data": Pxx, "centering": ["primal", "primal"]}, - {"name": popname + "_Pxy", "data": Pxy, "centering": ["primal", "primal"]}, - {"name": popname + "_Pxz", "data": Pxz, "centering": ["primal", "primal"]}, - {"name": popname + "_Pyy", "data": Pyy, "centering": ["primal", "primal"]}, - {"name": popname + "_Pyz", "data": Pyz, "centering": ["primal", "primal"]}, - {"name": popname + "_Pzz", "data": Pzz, "centering": ["primal", "primal"]}, + {"name": popname + "_Pxx", "data": Mxx.copy_as(Pxx)}, + {"name": popname + "_Pxy", "data": Mxy.copy_as(Pxy)}, + {"name": popname + "_Pxz", "data": Mxz.copy_as(Pxz)}, + {"name": popname + "_Pyy", "data": Myy.copy_as(Pyy)}, + {"name": popname + "_Pyz", "data": Myz.copy_as(Pyz)}, + {"name": popname + "_Pzz", "data": Mzz.copy_as(Pzz)}, ) diff --git a/pyphare/pyphare_tests/test_pharesee/test_geometry.py b/pyphare/pyphare_tests/test_pharesee/test_geometry.py index 6023f429f..c9ae8d1f5 100644 --- a/pyphare/pyphare_tests/test_pharesee/test_geometry.py +++ b/pyphare/pyphare_tests/test_pharesee/test_geometry.py @@ -116,7 +116,6 @@ def test_particle_ghost_area_boxes(self): } gaboxes = ghost_area_boxes(hierarchy, "particles") - particles = "particles" # same number of levels self.assertEqual(len(expected), len(gaboxes)) diff --git a/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py index 908484607..edd8a7fbd 100644 --- a/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py +++ b/pyphare/pyphare_tests/test_pharesee/test_hierarchy.py @@ -2,12 +2,12 @@ from ddt import ddt import numpy as np -from pyphare.simulator.simulator import Simulator +from pyphare.core.box import Box from pyphare.pharesee.run import Run +from pyphare.simulator.simulator import Simulator 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 @@ -152,12 +152,12 @@ def vthz(x, y): Simulator(config()).run() - def test_data_is_a_hierarchy(self): + def _test_data_is_a_hierarchy(self): r = Run(self.diag_dir()) B = r.GetB(0.0) self.assertTrue(isinstance(B, PatchHierarchy)) - def test_can_read_multiple_times(self): + def _test_can_read_multiple_times(self): r = Run(self.diag_dir()) times = (0.0, 0.1) B = r.GetB(times) @@ -168,7 +168,7 @@ def test_can_read_multiple_times(self): self.assertEqual(len(hier.times()), 2) self.assertTrue(np.allclose(hier.times().astype(np.float32), times)) - def test_hierarchy_is_refined(self): + def _test_hierarchy_is_refined(self): r = Run(self.diag_dir()) time = 0.0 B = r.GetB(time) @@ -177,19 +177,19 @@ def test_hierarchy_is_refined(self): self.assertEqual(len(B.levels(time)), 2) self.assertEqual(len(B.levels(time)), B.levelNbr(time)) - def test_can_get_nbytes(self): + def _test_can_get_nbytes(self): r = Run(self.diag_dir()) time = 0.0 B = r.GetB(time) self.assertGreater(B.nbytes(), 0) - def test_hierarchy_has_patches(self): + def _test_hierarchy_has_patches(self): r = Run(self.diag_dir()) time = 0.0 B = r.GetB(time) self.assertGreater(B.nbrPatches(), 0) - def test_access_patchdatas_as_hierarchies(self): + def _test_access_patchdatas_as_hierarchies(self): r = Run(self.diag_dir()) time = 0.0 B = r.GetB(time) @@ -197,7 +197,7 @@ def test_access_patchdatas_as_hierarchies(self): self.assertTrue(isinstance(B.y, PatchHierarchy)) self.assertTrue(isinstance(B.z, PatchHierarchy)) - def test_partial_domain_hierarchies(self): + def _test_partial_domain_hierarchies(self): import matplotlib.pyplot as plt from matplotlib.patches import Rectangle @@ -226,7 +226,7 @@ def test_partial_domain_hierarchies(self): self.assertTrue(isinstance(Bpartial, PatchHierarchy)) self.assertLess(Bpartial.nbrPatches(), B.nbrPatches()) - def test_scalarfield_quantities(self): + def _test_scalarfield_quantities(self): r = Run(self.diag_dir()) time = 0.0 Ni = r.GetNi(time) @@ -234,7 +234,7 @@ def test_scalarfield_quantities(self): self.assertTrue(isinstance(Ni, ScalarField)) self.assertTrue(isinstance(Vi, VectorField)) - def test_sum_two_scalarfields(self): + def _test_sum_two_scalarfields(self): r = Run(self.diag_dir()) time = 0.0 Ni = r.GetNi(time) @@ -243,7 +243,7 @@ def test_sum_two_scalarfields(self): self.assertTrue(isinstance(s, ScalarField)) self.assertEqual(s.quantities(), ["value"]) - def test_sum_with_scalar(self): + def _test_sum_with_scalar(self): r = Run(self.diag_dir()) time = 0.0 Ni = r.GetNi(time) @@ -253,7 +253,7 @@ def test_sum_with_scalar(self): self.assertTrue(isinstance(s, ScalarField)) self.assertEqual(s.quantities(), ["value"]) - def test_scalarfield_difference(self): + def _test_scalarfield_difference(self): r = Run(self.diag_dir()) time = 0.0 Ni = r.GetNi(time) @@ -265,7 +265,7 @@ def test_scalarfield_difference(self): self.assertTrue(isinstance(s, ScalarField)) self.assertEqual(s.quantities(), ["value"]) - def test_scalarfield_product(self): + def _test_scalarfield_product(self): r = Run(self.diag_dir()) time = 0.0 Ni = r.GetNi(time) @@ -276,7 +276,7 @@ def test_scalarfield_product(self): self.assertTrue(isinstance(s, ScalarField)) self.assertEqual(s.quantities(), ["value"]) - def test_scalarfield_division(self): + def _test_scalarfield_division(self): r = Run(self.diag_dir()) time = 0.0 Ni = r.GetNi(time) @@ -289,14 +289,14 @@ def test_scalarfield_division(self): self.assertTrue(isinstance(s, ScalarField)) self.assertEqual(s.quantities(), ["value"]) - def test_scalarfield_sqrt(self): + def _test_scalarfield_sqrt(self): r = Run(self.diag_dir()) time = 0.0 Ni = r.GetNi(time) nisqrt = sqrt(Ni) self.assertTrue(isinstance(nisqrt, ScalarField)) - def test_vectorfield_dot_product(self): + def _test_vectorfield_dot_product(self): r = Run(self.diag_dir()) time = 0.0 B = r.GetB(time) @@ -306,7 +306,7 @@ def test_vectorfield_dot_product(self): self.assertTrue(isinstance(s, ScalarField)) self.assertEqual(s.quantities(), ["value"]) - def test_vectorfield_binary_ops(self): + def _test_vectorfield_binary_ops(self): r = Run(self.diag_dir()) time = 0.0 B = r.GetB(time) @@ -321,7 +321,7 @@ def test_vectorfield_binary_ops(self): self.assertTrue(isinstance(s, VectorField)) self.assertEqual(s.quantities(), ["x", "y", "z"]) - def test_scalarfield_to_vectorfield_ops(self): + def _test_scalarfield_to_vectorfield_ops(self): r = Run(self.diag_dir()) time = 0.0 Ni = r.GetNi(time) @@ -330,7 +330,7 @@ def test_scalarfield_to_vectorfield_ops(self): self.assertTrue(isinstance(s, VectorField)) self.assertEqual(s.quantities(), ["x", "y", "z"]) - def test_vectorfield_to_vectorfield_ops(self): + def _test_vectorfield_to_vectorfield_ops(self): r = Run(self.diag_dir()) time = 0.0 B = r.GetB(time) @@ -340,6 +340,17 @@ def test_vectorfield_to_vectorfield_ops(self): self.assertTrue(isinstance(s, VectorField)) self.assertEqual(s.quantities(), ["x", "y", "z"]) + def test_all(self): + """ + DO NOT RUN MULTIPLE SIMULATIONS! + """ + + checks = 0 + for test in [method for method in dir(self) if method.startswith("_test_")]: + getattr(self, test)() + checks += 1 + self.assertEqual(checks, 18) # update if you add new tests + if __name__ == "__main__": unittest.main() diff --git a/src/amr/amr_constants.hpp b/src/amr/amr_constants.hpp index a7c43570d..4f263d754 100644 --- a/src/amr/amr_constants.hpp +++ b/src/amr/amr_constants.hpp @@ -1,10 +1,16 @@ -#ifndef AMR_CONSTANTS_HPP -#define AMR_CONSTANTS_HPP +#ifndef PHARE_AMR_CONSTANTS_HPP +#define PHARE_AMR_CONSTANTS_HPP #include -namespace PHARE::amr { + +namespace PHARE::amr +{ + static std::size_t constexpr refinementRatio = 2; -} + +static std::size_t constexpr MAX_LEVEL = 10; + +} // namespace PHARE::amr -#endif // AMR_CONSTANTS_HPP +#endif // PHARE_AMR_CONSTANTS_HPP diff --git a/src/amr/physical_models/hybrid_model.hpp b/src/amr/physical_models/hybrid_model.hpp index ecdcf1c9a..affbbd2e9 100644 --- a/src/amr/physical_models/hybrid_model.hpp +++ b/src/amr/physical_models/hybrid_model.hpp @@ -28,15 +28,15 @@ class HybridModel : public IPhysicalModel using Interface = IPhysicalModel; using amr_types = AMR_Types; using electrons_t = Electrons; - using patch_t = typename AMR_Types::patch_t; - using level_t = typename AMR_Types::level_t; + using patch_t = AMR_Types::patch_t; + using level_t = AMR_Types::level_t; using gridlayout_type = GridLayoutT; using electromag_type = Electromag; - using vecfield_type = typename Electromag::vecfield_type; - using field_type = typename vecfield_type::field_type; + using vecfield_type = Electromag::vecfield_type; + using field_type = vecfield_type::field_type; using grid_type = Grid_t; using ions_type = Ions; - using particle_array_type = typename Ions::particle_array_type; + using particle_array_type = Ions::particle_array_type; using resources_manager_type = amr::ResourcesManager; using ParticleInitializerFactory = core::ParticleInitializerFactory; diff --git a/src/amr/resources_manager/amr_utils.hpp b/src/amr/resources_manager/amr_utils.hpp index 48eb36cce..5ef0a2c01 100644 --- a/src/amr/resources_manager/amr_utils.hpp +++ b/src/amr/resources_manager/amr_utils.hpp @@ -4,8 +4,10 @@ #include "core/def.hpp" #include "core/def/phare_mpi.hpp" // IWYU pragma: keep +#include "core/utilities/mpi_utils.hpp" #include "core/utilities/constants.hpp" +#include "amr/amr_constants.hpp" #include "amr/types/amr_types.hpp" #include "amr/utilities/box/amr_box.hpp" @@ -16,6 +18,7 @@ #include #include #include +#include namespace PHARE { @@ -248,10 +251,69 @@ namespace amr iLevel++) { visitLevel(*hierarchy.getPatchLevel(iLevel), resman, - std::forward(action), std::forward(args...)); + std::forward(action), std::forward(args)...); } } + + template + auto boxesPerRankOn(auto const& level) + { + auto const& mapping = level.getProcessorMapping(); + auto const& global_boxes = level.getBoxes(); + + std::vector>> boxes_per_rank(core::mpi::size()); + assert(global_boxes.size() == level.getGlobalNumberOfPatches()); + + auto gbox_iter = global_boxes.begin(); + for (int i = 0; i < level.getGlobalNumberOfPatches(); ++i, ++gbox_iter) + boxes_per_rank[mapping.getProcessorAssignment(i)].emplace_back( + phare_box_from(*gbox_iter)); + + return boxes_per_rank; + } + + template + auto boxesPerRankOn(auto const& hierarchy, int const ilvl) + { + return boxesPerRankOn(hierarchy.getPatchLevel(ilvl)); + } + + + template + void onLevels(auto& hierarchy, Action&& action, std::size_t const minlvl = 0, + std::size_t const maxlvl = MAX_LEVEL) + { + if (hierarchy.getNumberOfLevels() < 1) + throw std::runtime_error("Hierarchy must have a level"); + + std::size_t const hier_levels = hierarchy.getNumberOfLevels() - 1; // size vs index + std::size_t const max = hier_levels < maxlvl ? hier_levels : maxlvl; + + for (std::size_t ilvl = minlvl; ilvl <= max; ++ilvl) + if (auto lvl = hierarchy.getPatchLevel(ilvl)) + action(*lvl); + } + + + + void onLevels(auto& hierarchy, auto&& onLevel, auto&& orMissing, std::size_t const minlvl, + std::size_t const maxlvl) + { + if (hierarchy.getNumberOfLevels() < 1) + throw std::runtime_error("Hierarchy must have a level"); + + std::size_t const hier_levels = hierarchy.getNumberOfLevels(); + + for (std::size_t ilvl = minlvl; ilvl < hier_levels; ++ilvl) + onLevel(*hierarchy.getPatchLevel(ilvl)); + + for (std::size_t ilvl = std::max(minlvl, hier_levels); ilvl <= maxlvl; ++ilvl) + orMissing(ilvl); + } + + + } // namespace amr } // namespace PHARE diff --git a/src/amr/wrappers/hierarchy.hpp b/src/amr/wrappers/hierarchy.hpp index aec20254b..53f0ecf1f 100644 --- a/src/amr/wrappers/hierarchy.hpp +++ b/src/amr/wrappers/hierarchy.hpp @@ -3,6 +3,7 @@ #include +#include "amr/amr_constants.hpp" #include "core/def/phare_mpi.hpp" // IWYU pragma: keep #include @@ -94,6 +95,7 @@ class Hierarchy : public HierarchyRestarter, public SAMRAI::hier::PatchHierarchy NO_DISCARD auto const& boundaryConditions() const { return boundaryConditions_; } NO_DISCARD auto const& cellWidth() const { return cellWidth_; } NO_DISCARD auto const& domainBox() const { return domainBox_; } + NO_DISCARD auto const& maxLevel() const { return maxLevel_; } @@ -120,6 +122,7 @@ class Hierarchy : public HierarchyRestarter, public SAMRAI::hier::PatchHierarchy std::vector const cellWidth_; std::vector const domainBox_; std::vector boundaryConditions_; + std::size_t maxLevel_ = 0; }; @@ -210,7 +213,16 @@ Hierarchy::Hierarchy(initializer::PHAREDict const& dict, , cellWidth_(cellWidth.data(), cellWidth.data() + dimension) , domainBox_(domainBox.data(), domainBox.data() + dimension) , boundaryConditions_(boundaryConditions.data(), boundaryConditions.data() + dimension) + { + auto const max_nbr_levels = dict["simulation"]["AMR"]["max_nbr_levels"].template to(); + if (max_nbr_levels < 1) + throw std::runtime_error("Invalid max_nbr_levels, must be >= 1"); + + maxLevel_ = max_nbr_levels - 1; + + if (maxLevel_ >= MAX_LEVEL) + throw std::runtime_error("Invalid max_nbr_levels, must be <= " + std::to_string(MAX_LEVEL)); } diff --git a/src/core/data/grid/gridlayout.hpp b/src/core/data/grid/gridlayout.hpp index 0aee24142..655d1f38b 100644 --- a/src/core/data/grid/gridlayout.hpp +++ b/src/core/data/grid/gridlayout.hpp @@ -2,19 +2,17 @@ #define PHARE_CORE_GRID_GridLayout_HPP -#include "core/hybrid/hybrid_quantities.hpp" +#include "core/def.hpp" #include "core/utilities/types.hpp" -#include "core/data/field/field.hpp" -#include "gridlayoutdefs.hpp" -#include "core/utilities/algorithm.hpp" #include "core/utilities/box/box.hpp" #include "core/utilities/constants.hpp" #include "core/utilities/index/index.hpp" #include "core/utilities/point/point.hpp" -#include "core/def.hpp" +#include "core/hybrid/hybrid_quantities.hpp" + +#include "gridlayoutdefs.hpp" #include -#include #include #include #include @@ -879,7 +877,6 @@ namespace core return GridLayoutImpl::centering(hybridQuantity); } - NO_DISCARD constexpr static std::array, 6> centering(HybridQuantity::Tensor hybridQuantity) { @@ -1027,6 +1024,29 @@ namespace core + /** + * @brief BxToMoments return the indexes and associated coef to compute the linear + * interpolation necessary to project Bx onto moments. + */ + NO_DISCARD auto static constexpr BxToMoments() { return GridLayoutImpl::BxToMoments(); } + + + /** + * @brief ByToMoments return the indexes and associated coef to compute the linear + * interpolation necessary to project By onto moments. + */ + NO_DISCARD auto static constexpr ByToMoments() { return GridLayoutImpl::ByToMoments(); } + + + /** + * @brief BzToMoments return the indexes and associated coef to compute the linear + * interpolation necessary to project Bz onto moments. + */ + NO_DISCARD auto static constexpr BzToMoments() { return GridLayoutImpl::BzToMoments(); } + + + + /** * @brief ExToMoments return the indexes and associated coef to compute the linear * interpolation necessary to project Ex onto moments. @@ -1157,30 +1177,23 @@ namespace core - - auto AMRGhostBoxFor(auto const& field) const + auto AMRBoxFor(auto const& field) const { auto const centerings = centering(field); - auto const growBy = [&]() { - std::array arr; - for (std::uint8_t i = 0; i < dimension; ++i) - arr[i] = nbrGhosts(centerings[i]); - return arr; - }(); - auto ghostBox = grow(AMRBox_, growBy); + auto box = AMRBox_; for (std::uint8_t i = 0; i < dimension; ++i) - ghostBox.upper[i] += (centerings[i] == QtyCentering::primal) ? 1 : 0; - return ghostBox; + box.upper[i] += (centerings[i] == QtyCentering::primal) ? 1 : 0; + return box; } - auto AMRBoxFor(auto const& field) const + auto AMRGhostBoxFor(auto const& field) const { auto const centerings = centering(field); - return grow(AMRGhostBoxFor(field), for_N_make_array([&](auto i) { - return -1 * nbrGhosts(centerings[i]); - })); + return grow(AMRBoxFor(field), for_N_make_array( + [&](auto i) { return nbrGhosts(centerings[i]); })); } + template void evalOnBox(Field& field, Fn&& fn) const { @@ -1203,6 +1216,15 @@ namespace core auto levelNumber() const { return levelNumber_; } + + auto amr_lcl_idx(auto const& box) const { return boxes_iterator{box, AMRToLocal(box)}; } + auto amr_lcl_idx() const { return amr_lcl_idx(AMRBox()); } + auto domain_amr_lcl_idx(auto const& field) const { return amr_lcl_idx(AMRBoxFor(field)); } + auto ghost_amr_lcl_idx(auto const& field) const + { + return amr_lcl_idx(AMRGhostBoxFor(field)); + } + private: template static void evalOnBox_(Field& field, Fn& fn, IndicesFn& startToEnd) @@ -1515,51 +1537,6 @@ namespace core } - - struct AMRLocalIndexer // - { - GridLayout const* layout; - Box amr_box; - Box lcl_box = layout->AMRToLocal(amr_box); - - struct iterator - { - void operator++() - { - ++amr_it; - ++lcl_it; - } - auto operator*() { return std::forward_as_tuple(*amr_it, *lcl_it); } - auto operator==(auto const& that) const - { - return amr_it == that.amr_it and lcl_it == that.lcl_it; - } - auto operator!=(auto const& that) const - { - return amr_it != that.amr_it or lcl_it != that.lcl_it; - } - - AMRLocalIndexer const* amr_lcl_indexer; - box_iterator amr_it; - box_iterator lcl_it; - }; - - auto begin() const { return iterator{this, amr_box.begin(), lcl_box.begin()}; } - auto end() const { return iterator{this, amr_box.end(), lcl_box.end()}; } - }; - - public: - auto amr_lcl_idx(auto const& box) const { return AMRLocalIndexer{this, box}; } - auto amr_lcl_idx() const { return amr_lcl_idx(AMRBox()); } - auto domain_amr_lcl_idx(auto const& field) const - { - return AMRLocalIndexer{this, AMRBoxFor(field)}; - } - auto ghost_amr_lcl_idx(auto const& field) const - { - return AMRLocalIndexer{this, AMRGhostBoxFor(field)}; - } - private: std::array meshSize_; Point origin_; diff --git a/src/core/data/grid/gridlayoutdefs.hpp b/src/core/data/grid/gridlayoutdefs.hpp index 23920396e..1b3cb02a0 100644 --- a/src/core/data/grid/gridlayoutdefs.hpp +++ b/src/core/data/grid/gridlayoutdefs.hpp @@ -1,11 +1,11 @@ #ifndef PHARE_CORE_GRID_GRIDLAYOUTDEFS_HPP #define PHARE_CORE_GRID_GRIDLAYOUTDEFS_HPP -#include -#include "core/hybrid/hybrid_quantities.hpp" -#include "core/utilities/types.hpp" #include "core/utilities/point/point.hpp" +#include "core/hybrid/hybrid_quantities.hpp" + +#include namespace PHARE { diff --git a/src/core/data/grid/gridlayoutimplyee.hpp b/src/core/data/grid/gridlayoutimplyee.hpp index d60537315..9a3746fef 100644 --- a/src/core/data/grid/gridlayoutimplyee.hpp +++ b/src/core/data/grid/gridlayoutimplyee.hpp @@ -424,19 +424,19 @@ namespace core { constexpr WeightPoint P1{Point{0}, 0.5}; constexpr WeightPoint P2{Point{iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 0.5}; constexpr WeightPoint P2{Point{iShift, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{iShift, 0, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -457,19 +457,19 @@ namespace core // in 1D the moment is already on Ey so return 1 point with no shift // with coef 1. constexpr WeightPoint P1{Point{0}, 1.}; - return std::array, 1>{P1}; + return std::array{P1}; } else if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, iShift, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -490,19 +490,116 @@ namespace core // in 1D or 2D the moment is already on Ez so return 1 point with no shift // with coef 1. constexpr WeightPoint P1{Point{0}, 1.}; - return std::array, 1>{P1}; + return std::array{P1}; } else if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 1.}; - return std::array, 1>{P1}; + return std::array{P1}; } else if constexpr (dimension == 3) { // in 3D we need two points, the second with a primalToDual shift along Z constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, 0, iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; + } + } + + + + NO_DISCARD auto static constexpr BxToMoments() + { + // Bx is primal dual dual + // moments are primal primal primal + // operation is thus Pdd to Ppp + [[maybe_unused]] auto constexpr iShift = dualToPrimal(); + + if constexpr (dimension == 1) + { + constexpr WeightPoint P1{Point{0}, 1.}; + return std::array{P1}; + } + if constexpr (dimension == 2) + { + constexpr WeightPoint P1{Point{0, 0}, 0.5}; + constexpr WeightPoint P2{Point{0, iShift}, 0.5}; + return std::array{P1, P2}; + } + else if constexpr (dimension == 3) + { + constexpr WeightPoint P1{Point{0, 0, 0}, 0.25}; + constexpr WeightPoint P2{Point{0, iShift, 0}, 0.25}; + constexpr WeightPoint P3{Point{0, 0, iShift}, 0.25}; + constexpr WeightPoint P4{Point{0, iShift, iShift}, 0.25}; + return std::array{P1, P2, P3, P4}; + } + } + + + + NO_DISCARD auto static constexpr ByToMoments() + { + // By is dual primal dual + // moments are primal primal primal + // operation is thus Dpd to Ppp + + auto constexpr iShift = dualToPrimal(); + + if constexpr (dimension == 1) + { + constexpr WeightPoint P1{Point{0}, 0.5}; + constexpr WeightPoint P2{Point{iShift}, 0.5}; + return std::array{P1, P2}; + } + if constexpr (dimension == 2) + { + constexpr WeightPoint P1{Point{0, 0}, 0.5}; + constexpr WeightPoint P2{Point{iShift, 0}, 0.5}; + return std::array{P1, P2}; + } + else if constexpr (dimension == 3) + { + constexpr WeightPoint P1{Point{0, 0, 0}, 0.25}; + constexpr WeightPoint P2{Point{iShift, 0, 0}, 0.25}; + constexpr WeightPoint P3{Point{0, 0, iShift}, 0.25}; + constexpr WeightPoint P4{Point{iShift, 0, iShift}, 0.25}; + return std::array{P1, P2, P3, P4}; + } + } + + + + + NO_DISCARD auto static constexpr BzToMoments() + { + // Bz is dual dual primal + // moments are primal primal primal + // operation is thus Ddp to Ppp + + auto constexpr iShift = dualToPrimal(); + + if constexpr (dimension == 1) + { + constexpr WeightPoint P1{Point{0}, 0.5}; + constexpr WeightPoint P2{Point{iShift}, 0.5}; + return std::array{P1, P2}; + } + if constexpr (dimension == 2) + { + constexpr WeightPoint P1{Point{0, 0}, 0.25}; + constexpr WeightPoint P2{Point{iShift, 0}, 0.25}; + constexpr WeightPoint P3{Point{0, iShift}, 0.25}; + constexpr WeightPoint P4{Point{iShift, iShift}, 0.25}; + return std::array{P1, P2, P3, P4}; + } + else if constexpr (dimension == 3) + { + constexpr WeightPoint P1{Point{0, 0, 0}, 0.25}; + constexpr WeightPoint P2{Point{iShift, 0, 0}, 0.25}; + constexpr WeightPoint P3{Point{0, iShift, 0}, 0.25}; + constexpr WeightPoint P4{Point{iShift, iShift, 0}, 0.25}; + return std::array{P1, P2, P3, P4}; } } @@ -521,19 +618,19 @@ namespace core { constexpr WeightPoint P1{Point{0}, 0.5}; constexpr WeightPoint P2{Point{iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 0.5}; constexpr WeightPoint P2{Point{iShift, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{iShift, 0, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -551,19 +648,19 @@ namespace core if constexpr (dimension == 1) { constexpr WeightPoint P1{Point{0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, iShift, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -581,18 +678,18 @@ namespace core if constexpr (dimension == 1) { constexpr WeightPoint P1{Point{0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, 0, iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -610,19 +707,19 @@ namespace core { constexpr WeightPoint P1{Point{0}, 0.5}; constexpr WeightPoint P2{Point{iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 0.5}; constexpr WeightPoint P2{Point{iShift, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{iShift, 0, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -639,19 +736,19 @@ namespace core if constexpr (dimension == 1) { constexpr WeightPoint P1{Point{0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, iShift, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -669,18 +766,18 @@ namespace core if constexpr (dimension == 1) { constexpr WeightPoint P1{Point{0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, 0, iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -696,7 +793,7 @@ namespace core { constexpr WeightPoint P1{Point{0}, 0.5}; constexpr WeightPoint P2{Point{p2dShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 2) { @@ -705,7 +802,7 @@ namespace core constexpr WeightPoint P3{Point{p2dShift, 0}, 0.25}; constexpr WeightPoint P4{Point{p2dShift, d2pShift}, 0.25}; - return std::array, 4>{P1, P2, P3, P4}; + return std::array{P1, P2, P3, P4}; } else if constexpr (dimension == 3) { @@ -722,13 +819,14 @@ namespace core 0.125}; constexpr WeightPoint P8{ Point{p2dShift, d2pShift, d2pShift}, 0.125}; - return std::array, 8>{P1, P2, P3, P4, P5, P6, P7, P8}; + return std::array{P1, P2, P3, P4, P5, P6, P7, P8}; } } NO_DISCARD auto static constexpr ByToEx() - { // By is dual primal dual + { + // By is dual primal dual // Ex is dual primal primal // operation is thus dpD to dpP // shift only in the Z direction @@ -737,18 +835,18 @@ namespace core if constexpr (dimension == 1) { constexpr WeightPoint P1{Point{0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, 0, iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -766,20 +864,20 @@ namespace core if constexpr (dimension == 1) { constexpr WeightPoint P1{Point{0}, 1.}; - return std::array, 1>{P1}; + return std::array{P1}; } else if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, iShift, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -796,7 +894,7 @@ namespace core { constexpr WeightPoint P1{Point{0}, 0.5}; constexpr WeightPoint P2{Point{d2pShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 2) @@ -806,7 +904,7 @@ namespace core constexpr WeightPoint P3{Point{0, d2pShift}, 0.25}; constexpr WeightPoint P4{Point{d2pShift, d2pShift}, 0.25}; - return std::array, 4>{P1, P2, P3, P4}; + return std::array{P1, P2, P3, P4}; } else if constexpr (dimension == 3) { @@ -822,7 +920,7 @@ namespace core 0.25}; constexpr WeightPoint P8{ Point{d2pShift, d2pShift, p2dShift}, 0.25}; - return std::array, 8>{P1, P2, P3, P4, P5, P6, P7, P8}; + return std::array{P1, P2, P3, P4, P5, P6, P7, P8}; } } @@ -841,19 +939,19 @@ namespace core { constexpr WeightPoint P1{Point{0}, 0.5}; constexpr WeightPoint P2{Point{iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 0.5}; constexpr WeightPoint P2{Point{iShift, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{iShift, 0, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -871,20 +969,20 @@ namespace core if constexpr (dimension == 1) { constexpr WeightPoint P1{Point{0}, 1.}; - return std::array, 1>{P1}; + return std::array{P1}; } else if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, iShift, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -902,19 +1000,19 @@ namespace core if constexpr (dimension == 1) { constexpr WeightPoint P1{Point{0}, 1}; - return std::array, 1>{P1}; + return std::array{P1}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 1.}; - return std::array, 1>{P1}; + return std::array{P1}; } else if constexpr (dimension == 3) { // in 3D we need two points, the second with a dualToPrimal shift along Z constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{0, 0, iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -931,7 +1029,7 @@ namespace core { constexpr WeightPoint P1{Point{0}, 0.5}; constexpr WeightPoint P2{Point{d2pShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 2) { @@ -940,7 +1038,7 @@ namespace core constexpr WeightPoint P3{Point{0, p2dShift}, 0.25}; constexpr WeightPoint P4{Point{d2pShift, p2dShift}, 0.25}; - return std::array, 4>{P1, P2, P3, P4}; + return std::array{P1, P2, P3, P4}; } else if constexpr (dimension == 3) { @@ -956,7 +1054,7 @@ namespace core 0.25}; constexpr WeightPoint P8{ Point{d2pShift, p2dShift, d2pShift}, 0.25}; - return std::array, 8>{P1, P2, P3, P4, P5, P6, P7, P8}; + return std::array{P1, P2, P3, P4, P5, P6, P7, P8}; } } @@ -973,19 +1071,19 @@ namespace core { constexpr WeightPoint P1{Point{0}, 0.5}; constexpr WeightPoint P2{Point{iShift}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 0.5}; constexpr WeightPoint P2{Point{iShift, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 0.5}; constexpr WeightPoint P2{Point{iShift, 0, 0}, 0.5}; - return std::array, 2>{P1, P2}; + return std::array{P1, P2}; } } @@ -1001,17 +1099,17 @@ namespace core if constexpr (dimension == 1) { constexpr WeightPoint P1{Point{0}, 1}; - return std::array, 1>{P1}; + return std::array{P1}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } } @@ -1027,17 +1125,17 @@ namespace core if constexpr (dimension == 1) { constexpr WeightPoint P1{Point{0}, 1}; - return std::array, 1>{P1}; + return std::array{P1}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } } @@ -1053,17 +1151,17 @@ namespace core if constexpr (dimension == 1) { constexpr WeightPoint P1{Point{0}, 1}; - return std::array, 1>{P1}; + return std::array{P1}; } if constexpr (dimension == 2) { constexpr WeightPoint P1{Point{0, 0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } else if constexpr (dimension == 3) { constexpr WeightPoint P1{Point{0, 0, 0}, 1.0}; - return std::array, 1>{P1}; + return std::array{P1}; } } }; // namespace core diff --git a/src/core/data/ndarray/ndarray_vector.hpp b/src/core/data/ndarray/ndarray_vector.hpp index a008c76ef..c4ab735a1 100644 --- a/src/core/data/ndarray/ndarray_vector.hpp +++ b/src/core/data/ndarray/ndarray_vector.hpp @@ -16,67 +16,63 @@ namespace PHARE::core { -template +template struct NdArrayViewer { - template - NO_DISCARD static DataType const& at(DataType const* data, NCells const& nCells, - Indexes const&... indexes) - { - auto params = std::forward_as_tuple(indexes...); - static_assert(sizeof...(Indexes) == dim); - // static_assert((... && std::is_unsigned_v)); TODO : manage later if - // this test should be included + using Id = std::uint32_t; + using Idx = Id const; + template typename Indexes, typename Index> + static inline std::uint32_t idx(NCells const& nCells, Indexes const& indexes) + { if constexpr (dim == 1) - { - auto i = std::get<0>(params); + return idx(nCells, indexes[0]); - return data[i]; - } + else if constexpr (dim == 2) + return idx(nCells, indexes[0], indexes[1]); - if constexpr (dim == 2) - { - auto i = std::get<0>(params); - auto j = std::get<1>(params); + else if constexpr (dim == 3) + return idx(nCells, indexes[0], indexes[1], indexes[2]); + } - if constexpr (c_ordering) - return data[j + i * nCells[1]]; - else - return data[i + j * nCells[0]]; - } + static inline std::uint32_t idx(auto const /*nCells*/, Idx i) { return i; } - if constexpr (dim == 3) - { - auto i = std::get<0>(params); - auto j = std::get<1>(params); - auto k = std::get<2>(params); - - if constexpr (c_ordering) - return data[k + j * nCells[2] + i * nCells[1] * nCells[2]]; - else - return data[i + j * nCells[0] + k * nCells[1] * nCells[0]]; - } + static inline std::uint32_t idx(auto const nCells, Idx i, Idx j) + { + if constexpr (c_ordering) + return j + i * nCells[1]; + else + return i + j * nCells[0]; + } + static inline std::uint32_t idx(auto const nCells, Idx i, Idx j, Idx k) + { + if constexpr (c_ordering) + return k + j * nCells[2] + i * nCells[1] * nCells[2]; + else + return i + j * nCells[0] + k * nCells[1] * nCells[0]; } - template typename Indexes, typename Index> - NO_DISCARD static DataType const& at(DataType const* data, NCells const& nCells, - Indexes const& indexes) + template typename Indexes, typename Index> + NO_DISCARD static inline auto& at(auto* data, auto const& nCells, + Indexes const& indexes) { - if constexpr (dim == 1) - return at(data, nCells, indexes[0]); - - else if constexpr (dim == 2) - return at(data, nCells, indexes[0], indexes[1]); + auto const& i = idx(nCells, indexes); + assert(i < product(nCells, std::uint32_t{1})); + return data[i]; + } - else if constexpr (dim == 3) - return at(data, nCells, indexes[0], indexes[1], indexes[2]); + static inline auto& at(auto* data, auto const nCells, auto const... indexes) + { + auto const& i = idx(nCells, indexes...); + assert(i < product(nCells, std::uint32_t{1})); + return data[i]; } }; + template class MaskedView { @@ -102,7 +98,7 @@ class MaskedView template NO_DISCARD DataType const& operator()(Indexes... indexes) const { - return NdArrayViewer::at(array_.data(), shape_, indexes...); + return NdArrayViewer::at(array_.data(), shape_, indexes...); } template @@ -139,7 +135,7 @@ class NdArrayView static std::size_t const dimension = dim; using type = DataType; using pointer_type = DataType*; - using viewer = NdArrayViewer; + using viewer = NdArrayViewer; explicit NdArrayView(pointer_type ptr, std::array const& nCells) : ptr_{ptr} @@ -292,7 +288,7 @@ class NdArrayVector template NO_DISCARD DataType const& operator()(Indexes... indexes) const { - return NdArrayViewer::at(data_.data(), nCells_, indexes...); + return NdArrayViewer::at(data_.data(), nCells_, indexes...); } template @@ -304,7 +300,7 @@ class NdArrayVector template NO_DISCARD DataType const& operator()(std::array const& indexes) const { - return NdArrayViewer::at(data_.data(), nCells_, indexes); + return NdArrayViewer::at(data_.data(), nCells_, indexes); } template diff --git a/src/core/data/tensorfield/tensorfield.hpp b/src/core/data/tensorfield/tensorfield.hpp index 48a6ca2d8..0a8afd270 100644 --- a/src/core/data/tensorfield/tensorfield.hpp +++ b/src/core/data/tensorfield/tensorfield.hpp @@ -205,7 +205,7 @@ class TensorField NO_DISCARD auto& componentNames() const { return componentNames_; } NO_DISCARD auto& physicalQuantity() const { return qty_; } - NO_DISCARD auto constexpr static size() { return N; } + NO_DISCARD auto static constexpr size() { return N; } private: auto static _get_index_for(Component component) diff --git a/src/core/utilities/algorithm.hpp b/src/core/utilities/algorithm.hpp index ae1aba764..286205264 100644 --- a/src/core/utilities/algorithm.hpp +++ b/src/core/utilities/algorithm.hpp @@ -1,13 +1,20 @@ -#ifndef PHARE_ALGORITHM_HPP -#define PHARE_ALGORITHM_HPP +#ifndef PHARE_CORE_UTILITIES_ALGORITHM_HPP +#define PHARE_CORE_UTILITIES_ALGORITHM_HPP + #include "core/def.hpp" -#include "core/utilities/span.hpp" +#include "core/utilities/span.hpp" +#include "core/utilities/types.hpp" +#include "core/data/field/field.hpp" +#include "core/utilities/box/box.hpp" +#include "core/data/grid/gridlayoutdefs.hpp" +#include "core/data/ndarray/ndarray_vector.hpp" +#include "core/data/tensorfield/tensorfield.hpp" -#include #include +#include namespace PHARE { @@ -64,6 +71,83 @@ void average(Span const& f1, Span const& f2, Span& avg) av[i] = (d1[i] + d2[i]) * .5; } + +template +auto convert_to_primal( // + auto const& src, // + auto const& layout, // + auto const lix, // + PhysicalQuantity const qty // +) +{ + using PQ = PhysicalQuantity; + + if (qty == PQ::Bx) + return layout.project(src, lix, layout.BxToMoments()); + else if (qty == PQ::By) + return layout.project(src, lix, layout.ByToMoments()); + else if (qty == PQ::Bz) + return layout.project(src, lix, layout.BzToMoments()); + + else if (qty == PQ::Ex) + return layout.project(src, lix, layout.ExToMoments()); + else if (qty == PQ::Ey) + return layout.project(src, lix, layout.EyToMoments()); + else if (qty == PQ::Ez) + return layout.project(src, lix, layout.EzToMoments()); + + throw std::runtime_error("Quantity not supported for conversion to primal."); +} + +template +auto& convert_to_fortran_primal( // DOES NOT WORK ON GHOST BOX! + Field& dst, // + Field const& src, // + auto const& layout // +) +{ + bool static constexpr c_ordering = false; + + assert(all(layout.centering(dst), [](auto const c) { return c == QtyCentering::primal; })); + + auto const all_primal + = all(layout.centering(src), [](auto const c) { return c == QtyCentering::primal; }); + + auto const lcl_box = layout.AMRToLocal(layout.AMRBoxFor(dst)); + auto dst_view = make_array_view(dst.data(), *lcl_box.shape()); + + auto const lcl_zero_box = [&]() { + auto box = lcl_box; + box.lower -= lcl_box.lower; // 0 + box.upper -= lcl_box.lower; // shift + return box; + }(); + + if (all_primal) + for (auto const [lix, lix0] : boxes_iterator(lcl_box, lcl_zero_box)) + dst_view(lix0) = src(lix); + + else + for (auto const [lix, lix0] : boxes_iterator(lcl_box, lcl_zero_box)) + dst_view(lix0) = convert_to_primal(src, layout, lix, src.physicalQuantity()); + + return dst; +} + +template +auto& convert_to_fortran_primal( // + TensorField& dst, // + TensorField const& src, // + auto const& layout // +) +{ + for (std::size_t ci = 0; ci < src.size(); ++ci) + convert_to_fortran_primal(dst[ci], src[ci], layout); + return dst; +} + + + } // namespace PHARE::core -#endif +#endif // PHARE_CORE_UTILITIES_ALGORITHM_HPP diff --git a/src/core/utilities/box/box.hpp b/src/core/utilities/box/box.hpp index 6e55d5196..a3efcf74c 100644 --- a/src/core/utilities/box/box.hpp +++ b/src/core/utilities/box/box.hpp @@ -2,15 +2,16 @@ #define PHARE_CORE_UTILITIES_BOX_BOX_HPP +#include "core/def.hpp" #include "core/utilities/types.hpp" #include "core/utilities/point/point.hpp" #include "core/utilities/meta/meta_utilities.hpp" -#include "core/def.hpp" +#include #include -#include #include #include +#include namespace PHARE::core { @@ -26,6 +27,7 @@ template struct Box { using value_type = Type; + using iterator = box_iterator; static constexpr auto dimension = dim; @@ -93,11 +95,11 @@ struct Box return *this; } - NO_DISCARD auto shape() const { return upper - lower + 1; } - NO_DISCARD auto size() const { return core::product(shape()); } + NO_DISCARD auto shape(std::size_t const i) const { return upper[i] - lower[i] + 1; } + NO_DISCARD auto shape() const { return upper - lower + 1; } + NO_DISCARD auto size() const { return core::product(shape(), std::size_t{1}); } - using iterator = box_iterator; NO_DISCARD auto begin() { return iterator{this, lower}; } // // since the 1D scan of the multidimensional box is done assuming C ordering @@ -188,6 +190,7 @@ class box_iterator return *this; } + bool operator!=(box_iterator const& other) const { return box_ != other.box_ or index_ != other.index_; @@ -203,6 +206,58 @@ class box_iterator template Box(Point lower, Point upper) -> Box; +template +struct boxes_iterator +{ + auto constexpr static N = sizeof...(Boxes); + + boxes_iterator(Boxes const&... boxes) + : boxes{std::forward_as_tuple(boxes...)} + { + } + + struct iterator + { + using Tuple_t = std::tuple; + static_assert(N == std::tuple_size_v); + + iterator(std::tuple iterators) + : its{iterators} + { + } + + void operator++() + { + for_N([&](auto i) { ++std::get(its); }); + } + + auto operator*() + { + return for_N([&](auto i) { return *std::get(its); }); + } + + auto operator!=(iterator const& that) const + { + return for_N_any([&](auto i) { return std::get(its) != std::get(that.its); }); + } + + Tuple_t its; + }; + + + + auto begin() + { + return iterator{for_N([&](auto i) { return std::get(boxes).begin(); })}; + } + auto end() + { + return iterator{for_N([&](auto i) { return std::get(boxes).end(); })}; + } + + std::tuple boxes; +}; + /** this overload of isIn takes a Point and a Container of boxes @@ -268,6 +323,7 @@ NO_DISCARD bool isIn(Point const& point, Box const& box) + template Box grow(Box const& box, OType const& size) { diff --git a/src/core/utilities/logger/logger_defaults.hpp b/src/core/utilities/logger/logger_defaults.hpp index 5ae689122..7f313604e 100644 --- a/src/core/utilities/logger/logger_defaults.hpp +++ b/src/core/utilities/logger/logger_defaults.hpp @@ -10,8 +10,8 @@ #endif // PHARE_WITH_PHLOP #ifndef PHARE_SCOPE_TIMER -#define PHARE_SCOPE_TIMER(str) // nothing -#endif // PHARE_SCOPE_TIMER +#define PHARE_SCOPE_TIMER(str) +#endif // PHARE_SCOPE_TIMER #if PHARE_LOG_LEVEL >= 1 diff --git a/src/diagnostic/detail/h5typewriter.hpp b/src/diagnostic/detail/h5typewriter.hpp index 57d450eb0..a05aea449 100644 --- a/src/diagnostic/detail/h5typewriter.hpp +++ b/src/diagnostic/detail/h5typewriter.hpp @@ -20,7 +20,7 @@ template class H5TypeWriter : public PHARE::diagnostic::TypeWriter { public: - using Attributes = typename Writer::Attributes; + using Attributes = Writer::Attributes; static constexpr auto dimension = Writer::dimension; H5TypeWriter(Writer& h5Writer) diff --git a/src/diagnostic/detail/h5writer.hpp b/src/diagnostic/detail/h5writer.hpp index 9e636e0b4..e6675d0d3 100644 --- a/src/diagnostic/detail/h5writer.hpp +++ b/src/diagnostic/detail/h5writer.hpp @@ -1,17 +1,22 @@ #ifndef PHARE_DETAIL_DIAGNOSTIC_HIGHFIVE_HPP #define PHARE_DETAIL_DIAGNOSTIC_HIGHFIVE_HPP - -#include "core/data/vecfield/vecfield_component.hpp" -#include "core/utilities/mpi_utils.hpp" #include "core/utilities/types.hpp" -#include "core/utilities/meta/meta_utilities.hpp" +#include "core/utilities/mpi_utils.hpp" +#include "core/data/vecfield/vecfield_component.hpp" + +#include "initializer/data_provider.hpp" #include "hdf5/detail/h5/h5_file.hpp" -#include "diagnostic/detail/h5typewriter.hpp" -#include "diagnostic/diagnostic_manager.hpp" + #include "diagnostic/diagnostic_props.hpp" +#include "diagnostic/detail/h5typewriter.hpp" +#include "diagnostic/detail/types/info.hpp" +#include "diagnostic/detail/types/meta.hpp" +#include "diagnostic/detail/types/fluid.hpp" +#include "diagnostic/detail/types/particle.hpp" +#include "diagnostic/detail/types/electromag.hpp" #if !defined(PHARE_DIAG_DOUBLES) @@ -23,17 +28,6 @@ namespace PHARE::diagnostic::h5 { using namespace hdf5::h5; -template -class ElectromagDiagnosticWriter; -template -class FluidDiagnosticWriter; -template -class ParticlesDiagnosticWriter; -template -class MetaDiagnosticWriter; -template -class InfoDiagnosticWriter; - template @@ -57,9 +51,9 @@ class H5Writer template H5Writer(Hierarchy& hier, Model& model, std::string const hifivePath, HiFile::AccessMode _flags) - : flags{_flags} + : modelView_{hier, model} + , flags{_flags} , filePath_{hifivePath} - , modelView_{hier, model} { } @@ -168,7 +162,11 @@ class H5Writer auto& modelView() { return modelView_; } auto timestamp() const { return timestamp_; } - std::size_t minLevel = 0, maxLevel = 10; // TODO hard-coded to be parametrized somehow +private: + ModelView modelView_; + +public: + std::size_t minLevel = 0, maxLevel = modelView_.maxLevel(); HiFile::AccessMode flags; @@ -176,7 +174,6 @@ class H5Writer double timestamp_ = 0; std::string filePath_; std::string patchPath_; // is passed around as "virtual write()" has no parameters - ModelView modelView_; Attributes fileAttributes_; std::unordered_map file_flags; @@ -201,8 +198,8 @@ class H5Writer H5Writer(H5Writer const&) = delete; H5Writer(H5Writer&&) = delete; - H5Writer& operator&(H5Writer const&) = delete; - H5Writer& operator&(H5Writer&&) = delete; + H5Writer& operator=(H5Writer const&) = delete; + H5Writer& operator=(H5Writer&&) = delete; // State of this class is controlled via "dump()" diff --git a/src/diagnostic/detail/vtk_types/electromag.hpp b/src/diagnostic/detail/vtk_types/electromag.hpp new file mode 100644 index 000000000..6865ea5cd --- /dev/null +++ b/src/diagnostic/detail/vtk_types/electromag.hpp @@ -0,0 +1,101 @@ +#ifndef PHARE_DIAGNOSTIC_DETAIL_VTK_TYPES_ELECTROMAG_HPP +#define PHARE_DIAGNOSTIC_DETAIL_VTK_TYPES_ELECTROMAG_HPP + +#include "diagnostic/detail/vtkh5_type_writer.hpp" + +#include +#include +#include + +namespace PHARE::diagnostic::vtkh5 +{ + +template +class ElectromagDiagnosticWriter : public H5TypeWriter +{ + using Super = H5TypeWriter; + using VTKFileWriter = Super::VTKFileWriter; + using VTKFileInitializer = Super::VTKFileInitializer; + +public: + ElectromagDiagnosticWriter(H5Writer& h5Writer) + : Super{h5Writer} + { + } + + void setup(DiagnosticProperties&) override; + void write(DiagnosticProperties&) override; + void compute(DiagnosticProperties&) override {} + +private: + struct Info + { + std::vector offset_per_level = std::vector(amr::MAX_LEVEL); + }; + + std::unordered_map mem; +}; + + +template +void ElectromagDiagnosticWriter::setup(DiagnosticProperties& diagnostic) +{ + auto& modelView = this->h5Writer_.modelView(); + VTKFileInitializer initializer{diagnostic, this}; + + if (mem.count(diagnostic.quantity) == 0) + mem.try_emplace(diagnostic.quantity); + auto& info = mem[diagnostic.quantity]; + + // assumes exists for all models + auto const init = [&](auto const ilvl) -> std::optional { + for (auto* vecField : this->h5Writer_.modelView().getElectromagFields()) + if (diagnostic.quantity == "/" + vecField->name()) + return initializer.template initTensorFieldFileLevel<1>(ilvl); + + return std::nullopt; + }; + + modelView.onLevels( + [&](auto const& level) { + auto const ilvl = level.getLevelNumber(); + if (auto const offset = init(ilvl)) + info.offset_per_level[ilvl] = *offset; + }, + [&](int const ilvl) { // missing level + init(ilvl); + }, + this->h5Writer_.minLevel, this->h5Writer_.maxLevel); +} + + + +template +void ElectromagDiagnosticWriter::write(DiagnosticProperties& diagnostic) +{ + auto& modelView = this->h5Writer_.modelView(); + auto& info = mem[diagnostic.quantity]; + + modelView.onLevels( + [&](auto const& level) { + auto const ilvl = level.getLevelNumber(); + + VTKFileWriter writer{diagnostic, this, info.offset_per_level[ilvl]}; + + auto const write_quantity = [&](auto& layout, auto const&, auto const) { + PHARE_LOG_SCOPE(3, "FluidDiagnosticWriter::write_quantity"); + + for (auto* vecField : this->h5Writer_.modelView().getElectromagFields()) + if (diagnostic.quantity == "/" + vecField->name()) + writer.template writeTensorField<1>(*vecField, layout); + }; + + modelView.visitHierarchy(write_quantity, ilvl, ilvl); + }, + this->h5Writer_.minLevel, this->h5Writer_.maxLevel); +} + + +} // namespace PHARE::diagnostic::vtkh5 + +#endif /* PHARE_DIAGNOSTIC_DETAIL_VTK_TYPES_ELECTROMAG_HPP */ diff --git a/src/diagnostic/detail/vtk_types/fluid.hpp b/src/diagnostic/detail/vtk_types/fluid.hpp new file mode 100644 index 000000000..4472716e4 --- /dev/null +++ b/src/diagnostic/detail/vtk_types/fluid.hpp @@ -0,0 +1,244 @@ +#ifndef PHARE_DIAGNOSTIC_DETAIL_VTK_TYPES_FLUID_HPP +#define PHARE_DIAGNOSTIC_DETAIL_VTK_TYPES_FLUID_HPP + +#include "core/logger.hpp" + +#include "amr/physical_models/mhd_model.hpp" +#include "amr/physical_models/hybrid_model.hpp" +#include "diagnostic/detail/vtkh5_type_writer.hpp" + +#include +#include +#include +#include + +namespace PHARE::diagnostic::vtkh5 +{ + +template +class FluidDiagnosticWriter : public H5TypeWriter +{ + using Super = H5TypeWriter; + using VTKFileWriter = Super::VTKFileWriter; + using VTKFileInitializer = Super::VTKFileInitializer; + +public: + FluidDiagnosticWriter(H5Writer& h5Writer) + : Super{h5Writer} + { + } + + void setup(DiagnosticProperties&) override; + void write(DiagnosticProperties&) override; + void compute(DiagnosticProperties&) override {}; + +private: + struct Info + { + std::vector offset_per_level = std::vector(amr::MAX_LEVEL); + }; + + struct HybridFluidInitializer + { + std::optional operator()(auto const ilvl); + + FluidDiagnosticWriter* writer; + DiagnosticProperties& diagnostic; + VTKFileInitializer& file_initializer; + }; + + struct MhdFluidInitializer + { + std::optional operator()(auto const ilvl); + + FluidDiagnosticWriter* writer; + DiagnosticProperties& diagnostic; + VTKFileInitializer& file_initializer; + }; + + struct HybridFluidWriter + { + void operator()(auto const& layout); + + FluidDiagnosticWriter* writer; + DiagnosticProperties& diagnostic; + VTKFileWriter& file_writer; + }; + + struct MhdFluidWriter + { + void operator()(auto const& layout); + + FluidDiagnosticWriter* writer; + DiagnosticProperties& diagnostic; + VTKFileWriter& file_writer; + }; + + auto static isActiveDiag(DiagnosticProperties const& diagnostic, std::string const& tree, + std::string const& var) + { + return diagnostic.quantity == tree + var; + }; + + std::unordered_map mem; +}; + + + +template +std::optional +FluidDiagnosticWriter::HybridFluidInitializer::operator()(auto const ilvl) +{ + auto& modelView = writer->h5Writer_.modelView(); + auto& ions = modelView.getIons(); + std::string const tree{"/ions/"}; + + if (isActiveDiag(diagnostic, tree, "charge_density")) + return file_initializer.initFieldFileLevel(ilvl); + if (isActiveDiag(diagnostic, tree, "mass_density")) + return file_initializer.initFieldFileLevel(ilvl); + if (isActiveDiag(diagnostic, tree, "bulkVelocity")) + return file_initializer.template initTensorFieldFileLevel<1>(ilvl); + + for (auto& pop : ions) + { + auto const pop_tree = tree + "pop/" + pop.name() + "/"; + if (isActiveDiag(diagnostic, pop_tree, "density")) + return file_initializer.initFieldFileLevel(ilvl); + else if (isActiveDiag(diagnostic, pop_tree, "charge_density")) + return file_initializer.initFieldFileLevel(ilvl); + else if (isActiveDiag(diagnostic, pop_tree, "flux")) + return file_initializer.template initTensorFieldFileLevel<1>(ilvl); + } + + return std::nullopt; +} + + + +template +std::optional +FluidDiagnosticWriter::MhdFluidInitializer::operator()(auto const ilvl) +{ + return std::nullopt; +} + + +template +void FluidDiagnosticWriter::setup(DiagnosticProperties& diagnostic) +{ + PHARE_LOG_SCOPE(3, "FluidDiagnosticWriter::setup"); + + using Model_t = H5Writer::ModelView::Model_t; + auto& modelView = this->h5Writer_.modelView(); + + VTKFileInitializer initializer{diagnostic, this}; + + if (mem.count(diagnostic.quantity) == 0) + mem.try_emplace(diagnostic.quantity); + auto& info = mem[diagnostic.quantity]; + + auto const init = [&](auto const ilvl) -> std::optional { + // + + if constexpr (solver::is_hybrid_model_v) + if (auto ret = HybridFluidInitializer{this, diagnostic, initializer}(ilvl)) + return ret; + + if constexpr (solver::is_mhd_model_v) + if (auto ret = MhdFluidInitializer{this, diagnostic, initializer}(ilvl)) + return ret; + + return std::nullopt; + }; + + modelView.onLevels( + [&](auto const& level) { + PHARE_LOG_SCOPE(3, "FluidDiagnosticWriter::setup_level"); + + auto const ilvl = level.getLevelNumber(); + if (auto const offset = init(ilvl)) + info.offset_per_level[ilvl] = *offset; + }, + [&](int const ilvl) { + PHARE_LOG_SCOPE(3, "FluidDiagnosticWriter::setup_missing_level"); + + init(ilvl); + }, + this->h5Writer_.minLevel, this->h5Writer_.maxLevel); +} + + + +template +void FluidDiagnosticWriter::HybridFluidWriter::operator()(auto const& layout) +{ + auto& modelView = writer->h5Writer_.modelView(); + auto& ions = modelView.getIons(); + std::string const tree{"/ions/"}; + + if (isActiveDiag(diagnostic, tree, "charge_density")) + file_writer.writeField(ions.chargeDensity(), layout); + else if (isActiveDiag(diagnostic, tree, "mass_density")) + file_writer.writeField(ions.massDensity(), layout); + else if (isActiveDiag(diagnostic, tree, "bulkVelocity")) + file_writer.template writeTensorField<1>(ions.velocity(), layout); + else + { + for (auto& pop : ions) + { + auto const pop_tree = tree + "pop/" + pop.name() + "/"; + if (isActiveDiag(diagnostic, pop_tree, "density")) + file_writer.writeField(pop.particleDensity(), layout); + else if (isActiveDiag(diagnostic, pop_tree, "charge_density")) + file_writer.writeField(pop.chargeDensity(), layout); + else if (isActiveDiag(diagnostic, pop_tree, "flux")) + file_writer.template writeTensorField<1>(pop.flux(), layout); + } + } +} + + + +template +void FluidDiagnosticWriter::MhdFluidWriter::operator()(auto const& layout) +{ + throw std::runtime_error("not implemented"); +} + + + +template +void FluidDiagnosticWriter::write(DiagnosticProperties& diagnostic) +{ + PHARE_LOG_SCOPE(3, "FluidDiagnosticWriter::write"); + + using Model_t = H5Writer::ModelView::Model_t; + auto& modelView = this->h5Writer_.modelView(); + auto const& info = mem[diagnostic.quantity]; + + modelView.onLevels( + [&](auto const& level) { + auto const ilvl = level.getLevelNumber(); + + VTKFileWriter writer{diagnostic, this, info.offset_per_level[ilvl]}; + + auto const write_quantity = [&](auto& layout, auto const&, auto const) { + if constexpr (solver::is_hybrid_model_v) + HybridFluidWriter{this, diagnostic, writer}(layout); + + if constexpr (solver::is_mhd_model_v) + MhdFluidWriter{this, diagnostic, writer}(layout); + }; + + modelView.visitHierarchy(write_quantity, ilvl, ilvl); + }, + this->h5Writer_.minLevel, this->h5Writer_.maxLevel); +} + + +} // namespace PHARE::diagnostic::vtkh5 + + + +#endif /* PHARE_DIAGNOSTIC_DETAIL_VTK_TYPES_FLUID_HPP */ diff --git a/src/diagnostic/detail/vtkh5_type_writer.hpp b/src/diagnostic/detail/vtkh5_type_writer.hpp new file mode 100644 index 000000000..9fd088350 --- /dev/null +++ b/src/diagnostic/detail/vtkh5_type_writer.hpp @@ -0,0 +1,466 @@ +#ifndef PHARE_DIAGNOSTIC_DETAIL_VTK_H5_TYPE_WRITER_HPP +#define PHARE_DIAGNOSTIC_DETAIL_VTK_H5_TYPE_WRITER_HPP + + +#include "core/logger.hpp" +#include "core/utilities/types.hpp" +#include "core/utilities/algorithm.hpp" +#include "core/utilities/mpi_utils.hpp" +#include "core/data/tensorfield/tensorfield.hpp" + +#include "amr/resources_manager/amr_utils.hpp" + +#include "diagnostic/diagnostic_writer.hpp" + +#include "hdf5/detail/h5/h5_file.hpp" + +#include +#include + +#if !defined(PHARE_DIAG_DOUBLES) +#error // PHARE_DIAG_DOUBLES not defined +#endif + + +namespace PHARE::diagnostic::vtkh5 +{ + +using namespace hdf5::h5; + + +template +struct HierarchyData +{ + // data is duplicated in lower dimensions to fit a 3d view + // so 2d data is * 2 + // and 1d data is * 4 + constexpr static auto X_TIMES = std::array{4, 2, 1}[dim - 1]; + + static auto& INSTANCE() + { + static HierarchyData data; + return data; + } + + static void reset(auto& h5Writer) // called once per dump + { + PHARE_LOG_SCOPE(3, "HierarchyData(level); + data.flattened_lcl_level_boxes[ilvl] + = flatten_boxes(data.level_boxes_per_rank[ilvl][core::mpi::rank()]); + + for (std::size_t i = 0; i < data.level_boxes_per_rank[ilvl].size(); ++i) + { + data.level_rank_data_size[ilvl].resize(core::mpi::size()); + + data.level_rank_data_size[ilvl][i] = 0; + for (auto box : data.level_boxes_per_rank[ilvl][i]) + { + box.upper += 1; // all primal + data.level_rank_data_size[ilvl][i] += box.size() * X_TIMES; // no ghosts + } + data.n_boxes_per_level[ilvl] += data.level_boxes_per_rank[ilvl][i].size(); + } + + data.level_data_size[ilvl] = core::sum(data.level_rank_data_size[ilvl]); + }, + [&](int const ilvl) { // ilvl does not exist currently + data.level_data_size[ilvl] = 0; + data.n_boxes_per_level[ilvl] = 0; + data.level_rank_data_size[ilvl].clear(); + data.level_rank_data_size[ilvl].resize(core::mpi::size()); + data.level_boxes_per_rank[ilvl].clear(); + data.level_boxes_per_rank[ilvl].resize(core::mpi::size()); + data.flattened_lcl_level_boxes[ilvl].clear(); + }, + h5Writer.minLevel, h5Writer.maxLevel); + } + + + auto static flatten_boxes(auto const& boxes) + { + std::vector> data(boxes.size()); + std::size_t pos = 0; + for (auto const& box : boxes) + { + for (std::size_t i = 0; i < dim; ++i) + { + data[pos][0 + i * 2] = box.lower[i]; + data[pos][1 + i * 2] = box.upper[i]; + } + ++pos; + } + return data; + } + + std::vector>> flattened_lcl_level_boxes; + std::vector>>> level_boxes_per_rank; + std::vector> level_rank_data_size; + std::vector n_boxes_per_level, level_data_size; +}; + + +template +class H5TypeWriter : public PHARE::diagnostic::TypeWriter +{ + using FloatType = std::conditional_t; + using HierData = HierarchyData; + using ModelView = Writer::ModelView; + using physical_quantity_type = ModelView::Field::physical_quantity_type; + using Attributes = Writer::Attributes; + std::string static inline const base = "/VTKHDF/"; + std::string static inline const level_base = base + "Level"; + std::string static inline const step_level = base + "Steps/Level"; + auto static inline const level_data_path + = [](auto const ilvl) { return level_base + std::to_string(ilvl) + "/PointData/data"; }; + + +public: + static constexpr auto dimension = Writer::dimension; + + H5TypeWriter(Writer& h5Writer) + : h5Writer_{h5Writer} + { + } + + virtual void setup(DiagnosticProperties&) = 0; + + void writeFileAttributes(DiagnosticProperties const&, Attributes const&); + +protected: + class VTKFileInitializer; + + class VTKFileWriter; + + + auto& getOrCreateH5File(DiagnosticProperties const& diagnostic) + { + if (!fileExistsFor(diagnostic)) + fileData_.emplace(diagnostic.quantity, this->h5Writer_.makeFile(diagnostic)); + return *fileData_.at(diagnostic.quantity); + } + + bool fileExistsFor(DiagnosticProperties const& diagnostic) const + { + return fileData_.count(diagnostic.quantity); + } + + + Writer& h5Writer_; + std::unordered_map> fileData_; +}; + + +template +void H5TypeWriter::writeFileAttributes(DiagnosticProperties const& prop, + Attributes const& attrs) +{ + if (fileExistsFor(prop)) // otherwise qty not supported (yet) + h5Writer_.writeGlobalAttributeDict(*fileData_.at(prop.quantity), attrs, "/"); +} + + +template +class H5TypeWriter::VTKFileInitializer +{ + auto static constexpr primal_qty = physical_quantity_type::rho; + std::size_t static constexpr boxValsIn3D = 6; // lo0, up0, lo1, up1, lo2, up2 +public: + VTKFileInitializer(DiagnosticProperties const& prop, H5TypeWriter* const typewriter); + + void initFileLevel(int const ilvl); + + std::size_t initFieldFileLevel(auto const ilvl) { return initAnyFieldLevel(ilvl); } + + template + std::size_t initTensorFieldFileLevel(auto const ilvl) + { + return initAnyFieldLevel()>(ilvl); + } + +private: + auto level_spacing(int const lvl) const + { + auto const mesh_size = typewriter->h5Writer_.modelView().cellWidth(); + return core::for_N_make_array( + [&](auto i) { return static_cast(mesh_size[i] / std::pow(2, lvl)); }); + } + + template + std::size_t initAnyFieldLevel(auto const ilvl) + { + h5file.create_resizable_2d_data_set(level_data_path(ilvl)); + initFileLevel(ilvl); + resize_boxes(ilvl); + resize_data(ilvl); + return data_offset; + } + + + template + void resize_data(int const ilvl); + + void resize_boxes(int const ilvl); + + bool const newFile; + DiagnosticProperties const& diagnostic; + H5TypeWriter* const typewriter; + HighFiveFile& h5file; + std::size_t data_offset; + ModelView& modelView = typewriter->h5Writer_.modelView(); +}; + + +template +class H5TypeWriter::VTKFileWriter +{ + auto static constexpr primal_qty = physical_quantity_type::rho; + std::size_t static constexpr boxValsIn3D = 6; // lo0, up0, lo1, up1, lo2, up2 + + +public: + VTKFileWriter(DiagnosticProperties const& prop, H5TypeWriter* const typewriter, + std::size_t const offset); + + void writeField(auto const& field, auto const& layout); + + template + void writeTensorField(auto const& tf, auto const& layout); + + + +private: + auto local_box(auto const& layout) const + { + return layout.AMRToLocal(layout.AMRBoxFor(primal_qty)); + } + + DiagnosticProperties const& diagnostic; + H5TypeWriter* const typewriter; + HighFiveFile& h5file; + std::size_t data_offset = 0; + ModelView& modelView = typewriter->h5Writer_.modelView(); +}; + + +template +H5TypeWriter::VTKFileInitializer::VTKFileInitializer(DiagnosticProperties const& prop, + H5TypeWriter* const tw) + : newFile{!tw->fileExistsFor(prop)} + , diagnostic{prop} + , typewriter{tw} + , h5file{tw->getOrCreateH5File(prop)} +{ + { + PHARE_LOG_SCOPE(3, "VTKFileInitializer::VTKFileInitializer::0"); + h5file.create_resizable_1d_data_set(base + "/Steps/Values"); + auto steps_group = h5file.file().getGroup(base + "/Steps"); + if (!steps_group.hasAttribute("NSteps")) + steps_group.template createAttribute("NSteps", HighFive::DataSpace::From(0)) + .write(0); + auto steps_attr = steps_group.getAttribute("NSteps"); + steps_attr.write(steps_attr.template read() + 1); + } + + { + PHARE_LOG_SCOPE(3, "VTKFileInitializer::VTKFileInitializer::1"); + auto const& timestamp = typewriter->h5Writer_.timestamp(); + auto ds = h5file.getDataSet(base + "/Steps/Values"); + auto const old_size = ds.getDimensions()[0]; + ds.resize({old_size + 1}); + ds.select({old_size}, {1}).write(timestamp); + } + + if (newFile) + { + PHARE_LOG_SCOPE(3, "VTKFileInitializer::VTKFileInitializer::2"); + auto root = h5file.file().getGroup(base); + + if (!root.hasAttribute("Version")) + root.template createAttribute>("Version", std::array{2, 2}); + + if (!root.hasAttribute("Origin")) + root.template createAttribute>("Origin", + std::array{0, 0, 0}); + + if (!root.hasAttribute("Type")) + root.template createAttribute("Type", "OverlappingAMR"); + + if (!root.hasAttribute("GridDescription")) + root.template createAttribute("GridDescription", "XYZ"); + } +} + + +template +H5TypeWriter::VTKFileWriter::VTKFileWriter(DiagnosticProperties const& prop, + H5TypeWriter* const tw, + std::size_t const offset) + : diagnostic{prop} + , typewriter{tw} + , h5file{tw->getOrCreateH5File(prop)} + , data_offset{offset} +{ +} + + +template +void H5TypeWriter::VTKFileWriter::writeField(auto const& field, auto const& layout) +{ + PHARE_LOG_SCOPE(3, "VTKFileWriter::writeField"); + + auto const& frimal = core::convert_to_fortran_primal(modelView.tmpField(), field, layout); + auto const size = local_box(layout).size(); + auto ds = h5file.getDataSet(level_data_path(layout.levelNumber())); + + PHARE_LOG_SCOPE(3, "VTKFileWriter::writeField::0"); + for (std::uint16_t i = 0; i < HierData::X_TIMES; ++i) + { + ds.select({data_offset, 0}, {size, 1}).write_raw(frimal.data()); + data_offset += size; + } +} + + +template +template +void H5TypeWriter::VTKFileWriter::writeTensorField(auto const& tf, auto const& layout) +{ + PHARE_LOG_SCOPE(3, "VTKFileWriter::writeTensorField"); + + auto static constexpr N = core::detail::tensor_field_dim_from_rank(); + auto const& frimal + = core::convert_to_fortran_primal(modelView.template tmpTensorField(), tf, layout); + auto const size = local_box(layout).size(); + auto ds = h5file.getDataSet(level_data_path(layout.levelNumber())); + + PHARE_LOG_SCOPE(3, "VTKFileWriter::writeTensorField::0"); + for (std::uint16_t i = 0; i < HierData::X_TIMES; ++i) + { + for (std::uint32_t c = 0; c < N; ++c) + ds.select({data_offset, c}, {size, 1}).write_raw(frimal[c].data()); + data_offset += size; + } +} + + +template +void H5TypeWriter::VTKFileInitializer::initFileLevel(int const ilvl) +{ + PHARE_LOG_SCOPE(3, "VTKFileInitializer::initFileLevel"); + + auto const lvl = std::to_string(ilvl); + + h5file.create_resizable_2d_data_set(level_base + lvl + "/AMRBox"); + h5file.create_resizable_1d_data_set(step_level + lvl + "/AMRBoxOffset"); + h5file.create_resizable_1d_data_set(step_level + lvl + "/NumberOfAMRBox"); + h5file.create_resizable_1d_data_set(step_level + lvl + "/PointDataOffset/data"); + + auto level_group = h5file.file().getGroup(level_base + lvl); + if (!level_group.hasAttribute("Spacing")) + { + PHARE_LOG_SCOPE(3, "VTKFileInitializer::initFileLevel::0"); + + level_group.template createAttribute>("Spacing", + level_spacing(ilvl)); + + level_group.createGroup("CellData"); + level_group.createGroup("FieldData"); + + auto steps_group = h5file.file().getGroup(step_level + lvl); + steps_group.createGroup("CellDataOffset"); + steps_group.createGroup("FieldDataOffset"); + } +} + + +template +template +void H5TypeWriter::VTKFileInitializer::resize_data(int const ilvl) +{ + PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_data"); + + auto const& hier_data = HierData::INSTANCE(); + auto const lvl = std::to_string(ilvl); + auto const data_path = level_data_path(ilvl); + auto point_data_ds = h5file.getDataSet(data_path); + data_offset = point_data_ds.getDimensions()[0]; + + { + PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_data::0"); + auto ds = h5file.getDataSet(step_level + lvl + "/PointDataOffset/data"); + auto const old_size = ds.getDimensions()[0]; + ds.resize({old_size + 1}); + ds.select({old_size}, {1}).write(data_offset); + } + + PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_data::1"); + auto const& rank_data_sizes = hier_data.level_rank_data_size[ilvl]; + auto const& level_data_size = hier_data.level_data_size[ilvl]; + auto const new_size = data_offset + level_data_size; + + point_data_ds.resize({new_size, N}); + for (int i = 0; i < core::mpi::rank(); ++i) + data_offset += rank_data_sizes[i]; +} + + +template +void H5TypeWriter::VTKFileInitializer::resize_boxes(int const ilvl) +{ + PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_boxes"); + + auto const lvl = std::to_string(ilvl); + auto const& hier_data = HierData::INSTANCE(); + auto const& rank_boxes = hier_data.level_boxes_per_rank[ilvl]; + auto const& total_boxes = hier_data.n_boxes_per_level[ilvl]; + + { + PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_boxes::0"); + auto ds = h5file.getDataSet(step_level + lvl + "/NumberOfAMRBox"); + auto const old_size = ds.getDimensions()[0]; + ds.resize({old_size + 1}); + ds.select({old_size}, {1}).write(total_boxes); + } + + auto amrbox_ds = h5file.getDataSet(level_base + lvl + "/AMRBox"); + auto box_offset = amrbox_ds.getDimensions()[0]; + + { + PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_boxes::1"); + auto ds = h5file.getDataSet(step_level + lvl + "/AMRBoxOffset"); + auto const old_size = ds.getDimensions()[0]; + ds.resize({old_size + 1}); + ds.select({old_size}, {1}).write(box_offset); + } + + PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_boxes::2"); + amrbox_ds.resize({box_offset + total_boxes, boxValsIn3D}); + for (int i = 0; i < core::mpi::rank(); ++i) + box_offset += rank_boxes[i].size(); + + + PHARE_LOG_SCOPE(3, "VTKFileInitializer::resize_boxes::3"); + amrbox_ds.select({box_offset, 0}, {rank_boxes[core::mpi::rank()].size(), dimension * 2}) + .write(hier_data.flattened_lcl_level_boxes[ilvl]); +} + + +} // namespace PHARE::diagnostic::vtkh5 + +#endif // PHARE_DIAGNOSTIC_DETAIL_VTK_H5_TYPE_WRITER_HPP diff --git a/src/diagnostic/detail/vtkh5_writer.hpp b/src/diagnostic/detail/vtkh5_writer.hpp new file mode 100644 index 000000000..e99b88b22 --- /dev/null +++ b/src/diagnostic/detail/vtkh5_writer.hpp @@ -0,0 +1,215 @@ +#ifndef PHARE_DIAGNOSTIC_DETAIL_VTK_H5_WRITER_HPP +#define PHARE_DIAGNOSTIC_DETAIL_VTK_H5_WRITER_HPP + + +#include "core/logger.hpp" +#include "core/utilities/mpi_utils.hpp" + +#include "amr/amr_constants.hpp" + +#include "initializer/data_provider.hpp" + +#include "diagnostic/diagnostic_props.hpp" +#include "diagnostic/detail/vtkh5_type_writer.hpp" + +#include "diagnostic/detail/vtk_types/fluid.hpp" +#include "diagnostic/detail/vtk_types/electromag.hpp" + +#include "hdf5/detail/h5/h5_file.hpp" + + +namespace PHARE::diagnostic::vtkh5 +{ +using namespace hdf5::h5; + + + +template +class H5Writer +{ + struct NullTypeWriter : public H5TypeWriter> + { + NullTypeWriter(auto& h5Writer) + : H5TypeWriter>{h5Writer} + { + } + + void setup(DiagnosticProperties& prop) {} + void write(DiagnosticProperties& prop) + { + if (core::mpi::rank() == 0) + { + PHARE_LOG_LINE_SS( // + "No diagnostic writer found for " + prop.type + ":" + prop.quantity); + } + } + void compute(DiagnosticProperties&) {} + }; + +public: + using ModelView = _ModelView; + using This = H5Writer; + using GridLayout = ModelView::GridLayout; + using Attributes = ModelView::PatchProperties; + + static constexpr auto dimension = GridLayout::dimension; + static constexpr auto interpOrder = GridLayout::interp_order; + static constexpr auto READ_WRITE = HiFile::AccessMode::OpenOrCreate; + + // flush_never: disables manual file closing, but still occurrs via RAII + static constexpr std::size_t flush_never = 0; + + template + H5Writer(Hierarchy& hier, Model& model, std::string const hifivePath, HiFile::AccessMode _flags) + : modelView_{hier, model} + , flags{_flags} + , filePath_{hifivePath} + { + } + + ~H5Writer() {} + H5Writer(H5Writer const&) = delete; + H5Writer(H5Writer&&) = delete; + H5Writer& operator=(H5Writer const&) = delete; + H5Writer& operator=(H5Writer&&) = delete; + + + template + static auto make_unique(Hierarchy& hier, Model& model, initializer::PHAREDict const& dict) + { + std::string filePath = dict["filePath"].template to(); + HiFile::AccessMode flags = READ_WRITE; + if (dict.contains("mode") and dict["mode"].template to() == "overwrite") + flags |= HiFile::Truncate; + return std::make_unique(hier, model, filePath, flags); + } + + + void dump(std::vector const&, double current_timestamp); + void dump_level(std::size_t level, std::vector const& diagnostics, + double timestamp); + + template + auto getDiagnosticWriterForType(String& type) + { + return typeWriters_.at(type); + } + + + auto makeFile(DiagnosticProperties const& diagnostic) + { + return std::make_unique(filePath_ + "/" + fileString(diagnostic.quantity), + file_flags[diagnostic.type + diagnostic.quantity]); + } + + template + static void writeGlobalAttributeDict(HighFiveFile& h5, Dict const& dict, std::string path) + { + dict.visit( + [&](std::string const& key, auto const& val) { h5.write_attribute(path, key, val); }); + } + + + auto& modelView() { return modelView_; } + auto timestamp() const { return timestamp_; } + +private: + ModelView modelView_; + +public: + std::size_t minLevel = 0, maxLevel = modelView_.maxLevel(); + HiFile::AccessMode flags; + + +private: + double timestamp_ = 0; + std::string filePath_; + Attributes fileAttributes_; + + std::unordered_map file_flags; + + std::unordered_map>> typeWriters_{ + {"info", make_writer()}, + {"meta", make_writer()}, + {"fluid", make_writer>()}, + {"electromag", make_writer>()}, + {"particle", make_writer()} // + }; + + template + std::shared_ptr> make_writer() + { + return std::make_shared(*this); + } + + + static std::string fileString(std::string fileStr) + { + if (fileStr[0] == '/') + fileStr = fileStr.substr(1); + std::replace(fileStr.begin(), fileStr.end(), '/', '_'); + return fileStr + ".vtkhdf"; + } + + + // State of this class is controlled via "dump()" + // block public access to internal state + friend class H5TypeWriter; + friend class FluidDiagnosticWriter; + friend class ElectromagDiagnosticWriter; +}; + + + +template +void H5Writer::dump(std::vector const& diagnostics, + double timestamp) +{ + timestamp_ = timestamp; + fileAttributes_["dimension"] = dimension; + fileAttributes_["interpOrder"] = interpOrder; + fileAttributes_["domain_box"] = modelView_.domainBox(); + fileAttributes_["boundary_conditions"] = modelView_.boundaryConditions(); + + HierarchyData::reset(*this); + + for (auto* diagnostic : diagnostics) + if (!file_flags.count(diagnostic->type + diagnostic->quantity)) + file_flags[diagnostic->type + diagnostic->quantity] = this->flags; + + for (auto* diagnostic : diagnostics) // all collective calls first! + { + auto& type_writer = *typeWriters_.at(diagnostic->type); + type_writer.setup(*diagnostic); + type_writer.writeFileAttributes(*diagnostic, fileAttributes_); + } + + for (auto* diagnostic : diagnostics) + typeWriters_.at(diagnostic->type)->write(*diagnostic); + + for (auto* diagnostic : diagnostics) // don't truncate past first dump + file_flags[diagnostic->type + diagnostic->quantity] = READ_WRITE; +} + +template +void H5Writer::dump_level(std::size_t level, + std::vector const& diagnostics, + double timestamp) +{ + std::size_t _minLevel = this->minLevel; + std::size_t _maxLevel = this->maxLevel; + + this->minLevel = level; + this->maxLevel = level; + + this->dump(diagnostics, timestamp); + + this->minLevel = _minLevel; + this->maxLevel = _maxLevel; +} + + + +} // namespace PHARE::diagnostic::vtkh5 + +#endif /* PHARE_DIAGNOSTIC_DETAIL_VTK_H5_WRITER_HPP */ diff --git a/src/diagnostic/diagnostic_manager.hpp b/src/diagnostic/diagnostic_manager.hpp index a7d0f0ea6..14c6560b2 100644 --- a/src/diagnostic/diagnostic_manager.hpp +++ b/src/diagnostic/diagnostic_manager.hpp @@ -31,7 +31,7 @@ void registerDiagnostics(DiagManager& dMan, initializer::PHAREDict const& diagsP while (diagsParams.contains(diagType) && diagsParams[diagType].contains(diagType + std::to_string(diagBlockID))) { - const std::string diagName = diagType + std::to_string(diagBlockID); + std::string const diagName = diagType + std::to_string(diagBlockID); dMan.addDiagDict(diagsParams[diagType][diagName]); ++diagBlockID; } diff --git a/src/diagnostic/diagnostic_model_view.hpp b/src/diagnostic/diagnostic_model_view.hpp index 4885db24d..91a1d11eb 100644 --- a/src/diagnostic/diagnostic_model_view.hpp +++ b/src/diagnostic/diagnostic_model_view.hpp @@ -1,6 +1,7 @@ #ifndef DIAGNOSTIC_MODEL_VIEW_HPP #define DIAGNOSTIC_MODEL_VIEW_HPP +#include "amr/amr_constants.hpp" #include "core/def.hpp" #include "core/utilities/mpi_utils.hpp" @@ -33,8 +34,8 @@ class BaseModelView : public IModelView using GridLayout = Model::gridlayout_type; using VecField = Model::vecfield_type; using TensorFieldT = Model::ions_type::tensorfield_type; - using GridLayoutT = Model::gridlayout_type; using ResMan = Model::resources_manager_type; + using Field = Model::field_type; using TensorFieldData_t = ResMan::template UserTensorField_t::patch_data_type; static constexpr auto dimension = Model::dimension; @@ -67,41 +68,50 @@ class BaseModelView : public IModelView auto& rm = *model_.resourcesManager; auto& ions = model_.state.ions; - for (auto patch : rm.enumerate(lvl, ions, sumTensor_)) + for (auto patch : rm.enumerate(lvl, ions, tmpTensor_)) for (std::uint8_t c = 0; c < N; ++c) - std::memcpy(sumTensor_[c].data(), ions[popidx].momentumTensor()[c].data(), + std::memcpy(tmpTensor_[c].data(), ions[popidx].momentumTensor()[c].data(), ions[popidx].momentumTensor()[c].size() * sizeof(value_type)); MTAlgos[popidx].getOrCreateSchedule(hierarchy_, lvl.getLevelNumber()).fillData(time); - for (auto patch : rm.enumerate(lvl, ions, sumTensor_)) + for (auto patch : rm.enumerate(lvl, ions, tmpTensor_)) for (std::uint8_t c = 0; c < N; ++c) - std::memcpy(ions[popidx].momentumTensor()[c].data(), sumTensor_[c].data(), + std::memcpy(ions[popidx].momentumTensor()[c].data(), tmpTensor_[c].data(), ions[popidx].momentumTensor()[c].size() * sizeof(value_type)); } template - void onLevels(Action&& action, int minlvl = 0, int maxlvl = 0) + void onLevels(Action&& action, std::size_t const minlvl = 0, + std::size_t const maxlvl = amr::MAX_LEVEL) { - for (int ilvl = minlvl; ilvl < hierarchy_.getNumberOfLevels() && ilvl <= maxlvl; ++ilvl) - if (auto lvl = hierarchy_.getPatchLevel(ilvl)) - action(*lvl); + amr::onLevels(hierarchy_, std::forward(action), minlvl, maxlvl); + } + + + template + void onLevels(OnLevel&& onLevel, OrMissing&& orMissing, std::size_t const minlvl, + std::size_t const maxlvl) + { + amr::onLevels(hierarchy_, std::forward(onLevel), + std::forward(orMissing), minlvl, maxlvl); } template void visitHierarchy(Action&& action, int minLevel = 0, int maxLevel = 0) { - PHARE::amr::visitHierarchy(hierarchy_, *model_.resourcesManager, - std::forward(action), minLevel, maxLevel, - model_); + amr::visitHierarchy(hierarchy_, *model_.resourcesManager, + std::forward(action), minLevel, maxLevel, *this, + model_); } NO_DISCARD auto boundaryConditions() const { return hierarchy_.boundaryConditions(); } NO_DISCARD auto domainBox() const { return hierarchy_.domainBox(); } NO_DISCARD auto origin() const { return std::vector(dimension, 0); } NO_DISCARD auto cellWidth() const { return hierarchy_.cellWidth(); } + NO_DISCARD auto maxLevel() const { return hierarchy_.maxLevel(); } NO_DISCARD std::string getLayoutTypeString() const { @@ -141,6 +151,37 @@ class BaseModelView : public IModelView return model_.tags.at(key); } + + // auto localLevelBoxes(auto const ilvl) + // { + // std::vector> boxes; + // auto const& lvl = *hierarchy_.getPatchLevel(ilvl); + // boxes.reserve(lvl.getLocalNumberOfPatches()); + // visitHierarchy( + // [&](auto& layout, auto const&, auto const) { boxes.emplace_back(layout.AMRBox()); }, + // ilvl, ilvl); + // return boxes; + // } + + auto& tmpField() { return tmpField_; } + auto& tmpVecField() { return tmpVec_; } + + template + auto& tmpTensorField() + { + static_assert(rank > 0 and rank < 3); + if constexpr (rank == 1) + return tmpVec_; + else + return tmpTensor_; + } + + NO_DISCARD auto getCompileTimeResourcesViewList() + { + return std::forward_as_tuple(tmpField_, tmpVec_, tmpTensor_); + } + + protected: Model& model_; Hierarchy& hierarchy_; @@ -149,7 +190,7 @@ class BaseModelView : public IModelView { auto& rm = *model_.resourcesManager; - auto const dst_name = sumTensor_.name(); + auto const dst_name = tmpTensor_.name(); for (auto& pop : model_.state.ions) { @@ -160,7 +201,7 @@ class BaseModelView : public IModelView MTAlgo.MTalgo->registerRefine( idDst, idSrc, idDst, nullptr, std::make_shared< - amr::TensorFieldGhostInterpOverlapFillPattern>()); + amr::TensorFieldGhostInterpOverlapFillPattern>()); } // can't create schedules here as the hierarchy has no levels yet @@ -185,7 +226,9 @@ class BaseModelView : public IModelView }; std::vector MTAlgos; - TensorFieldT sumTensor_{"PHARE_sumTensor", core::HybridQuantity::Tensor::M}; + Field tmpField_{"PHARE_sumField", core::HybridQuantity::Scalar::rho}; + VecField tmpVec_{"PHARE_sumVec", core::HybridQuantity::Vector::V}; + TensorFieldT tmpTensor_{"PHARE_sumTensor", core::HybridQuantity::Tensor::M}; }; diff --git a/src/diagnostic/diagnostic_writer.hpp b/src/diagnostic/diagnostic_writer.hpp index a33c0a403..1d0b22344 100644 --- a/src/diagnostic/diagnostic_writer.hpp +++ b/src/diagnostic/diagnostic_writer.hpp @@ -3,7 +3,7 @@ #include "diagnostic_props.hpp" -#include "dict.hpp" + namespace PHARE::diagnostic { diff --git a/src/diagnostic/diagnostics.hpp b/src/diagnostic/diagnostics.hpp index 74111ff37..041efc374 100644 --- a/src/diagnostic/diagnostics.hpp +++ b/src/diagnostic/diagnostics.hpp @@ -25,11 +25,7 @@ #include "diagnostic_model_view.hpp" #include "diagnostic/detail/h5writer.hpp" -#include "diagnostic/detail/types/electromag.hpp" -#include "diagnostic/detail/types/particle.hpp" -#include "diagnostic/detail/types/fluid.hpp" -#include "diagnostic/detail/types/meta.hpp" -#include "diagnostic/detail/types/info.hpp" +#include "diagnostic/detail/vtkh5_writer.hpp" #endif @@ -56,8 +52,13 @@ struct DiagnosticsManagerResolver { #if PHARE_HAS_HIGHFIVE using ModelView_t = ModelView; - using Writer_t = h5::H5Writer; - return DiagnosticsManager::make_unique(hier, model, dict); + auto const format = cppdict::get_value(dict, "format", std::string{"phareh5"}); + + if (format == "phareh5") + return DiagnosticsManager>::make_unique(hier, model, dict); + if (format == "pharevtkhdf") + return DiagnosticsManager>::make_unique(hier, model, dict); + throw std::runtime_error("DiagnosticsManagerResolver - unknown format " + format); #else return std::make_unique(); #endif diff --git a/src/hdf5/detail/h5/h5_file.hpp b/src/hdf5/detail/h5/h5_file.hpp index 510c0786b..f8208d75c 100644 --- a/src/hdf5/detail/h5/h5_file.hpp +++ b/src/hdf5/detail/h5/h5_file.hpp @@ -3,30 +3,29 @@ #include "core/def.hpp" #include "core/def/phare_mpi.hpp" // IWYU pragma: keep -#include "highfive/H5File.hpp" -#include "highfive/H5Easy.hpp" - #include "core/utilities/types.hpp" #include "core/utilities/mpi_utils.hpp" #include "core/utilities/meta/meta_utilities.hpp" -namespace PHARE::hdf5::h5 + +#include "highfive/H5File.hpp" +#include "highfive/H5Easy.hpp" + + + +namespace PHARE::hdf5::h5::detail { -using HiFile = HighFive::File; -using FileOp = HighFive::File::AccessMode; +// https://support.hdfgroup.org/documentation/hdf5/latest/hdf5_chunking.html +static inline auto const CHUNK_SIZE = core::get_env_as("PHARE_H5_CHUNK_SIZE", std::size_t{1024}); +} // namespace PHARE::hdf5::h5::detail -template -NO_DISCARD auto decay_to_pointer(Data& data) +namespace PHARE::hdf5::h5 { - if constexpr (dim == 1) - return data.data(); - if constexpr (dim == 2) - return data[0].data(); - if constexpr (dim == 3) - return data[0][0].data(); -} +using HiFile = HighFive::File; +using FileOp = HighFive::File::AccessMode; + template NO_DISCARD auto vector_for_dim() @@ -107,11 +106,41 @@ class HighFiveFile template auto create_data_set(std::string const& path, Size const& dataSetSize) { + if (exist(path)) + return h5file_.getDataSet(path); createGroupsToDataSet(path); return h5file_.createDataSet(path, HighFive::DataSpace(dataSetSize)); } + template + auto create_chunked_data_set(auto const& path, auto const chunk, auto const& dataspace) + { + if (exist(path)) + return h5file_.getDataSet(path); + createGroupsToDataSet(path); + HighFive::DataSetCreateProps props; + props.add(HighFive::Chunking{chunk}); + return h5file_.createDataSet(path, dataspace, HighFive::create_datatype(), props); + } + + template + auto create_resizable_2d_data_set(auto const& path, hsize_t const x_chunk = detail::CHUNK_SIZE, + hsize_t const y_chunk = 1) + { + return create_chunked_data_set( + path, std::vector{x_chunk, y_chunk}, + HighFive::DataSpace({0, cols}, {HighFive::DataSpace::UNLIMITED, cols})); + } + + template + auto create_resizable_1d_data_set(auto const& path, hsize_t const chunk = detail::CHUNK_SIZE) + { + return create_chunked_data_set( + path, std::vector{chunk}, + HighFive::DataSpace({0}, {HighFive::DataSpace::UNLIMITED})); + } + /* * Communicate all dataset paths and sizes to all MPI process to allow each to create all @@ -214,7 +243,7 @@ class HighFiveFile }(); auto const paths = core::mpi::collect(path, mpi_size); - for (int i = 0; i < mpi_size; i++) + for (int i = 0; i < mpi_size; ++i) { std::string const keyPath = paths[i] == "null" ? "" : paths[i]; if (keyPath.empty()) @@ -244,11 +273,18 @@ class HighFiveFile return attr; } + bool exist(std::string const& s) const { return h5file_.exist(s); } + auto getDataSet(std::string const& s) + { + if (!exist(s)) + throw std::runtime_error("Dataset does not exist: " + s); + return h5file_.getDataSet(s); + } - HighFiveFile(HighFiveFile const&) = delete; - HighFiveFile(HighFiveFile const&&) = delete; - HighFiveFile& operator=(HighFiveFile const&) = delete; - HighFiveFile& operator=(HighFiveFile const&&) = delete; + HighFiveFile(HighFiveFile const&) = delete; + HighFiveFile(HighFiveFile&&) = delete; + HighFiveFile& operator=(HighFiveFile const&) = delete; + HighFiveFile& operator=(HighFiveFile&&) = delete; private: HighFive::FileAccessProps fapl_; diff --git a/src/phare/CMakeLists.txt b/src/phare/CMakeLists.txt index 0dd7d5bd3..fe1a459f9 100644 --- a/src/phare/CMakeLists.txt +++ b/src/phare/CMakeLists.txt @@ -4,7 +4,6 @@ project(phare-exe) add_executable(phare-exe ${SOURCES_INC} phare.cpp) target_compile_options(phare-exe PRIVATE ${PHARE_WERROR_FLAGS} -DPHARE_HAS_HIGHFIVE=${PHARE_HAS_HIGHFIVE}) -configure_file(${CMAKE_CURRENT_SOURCE_DIR}/phare_init.py ${CMAKE_CURRENT_BINARY_DIR}/phare_init.py @ONLY) target_link_libraries(phare-exe PUBLIC ${PHARE_BASE_LIBS} phare_simulator diff --git a/src/phare/phare_init.py b/src/phare/phare_init.py deleted file mode 100644 index d79de42b7..000000000 --- a/src/phare/phare_init.py +++ /dev/null @@ -1,109 +0,0 @@ -#!/usr/bin/env python3 - -import numpy as np - -import pyphare.pharein as ph - -# configure the simulation - -ph.Simulation( - smallest_patch_size=20, - largest_patch_size=20, - time_step_nbr=30000, # number of time steps (not specified if time_step and final_time provided) - final_time=30.0, # simulation final time (not specified if time_step and time_step_nbr provided) - boundary_types="periodic", # boundary condition, string or tuple, length == len(cell) == len(dl) - cells=40, # integer or tuple length == dimension - dl=0.3, # mesh size of the root level, float or tuple - max_nbr_levels=2, # (default=1) max nbr of levels in the AMR hierarchy - refinement="tagging", - # refinement_boxes={"L0": {"B0": [(10, ), (20, )]}}, - diag_options={"format": "phareh5", "options": {"dir": "phare_outputs"}}, -) - - -def density(x): - return 1.0 - - -def by(x): - from pyphare.pharein.global_vars import sim - - L = sim.simulation_domain() - return 0.1 * np.cos(2 * np.pi * x / L[0]) - - -def bz(x): - from pyphare.pharein.global_vars import sim - - L = sim.simulation_domain() - return 0.1 * np.sin(2 * np.pi * x / L[0]) - - -def bx(x): - return 1.0 - - -def vx(x): - return 0.0 - - -def vy(x): - from pyphare.pharein.global_vars import sim - - L = sim.simulation_domain() - return 0.1 * np.cos(2 * np.pi * x / L[0]) - - -def vz(x): - from pyphare.pharein.global_vars import sim - - L = sim.simulation_domain() - return 0.1 * np.sin(2 * np.pi * x / L[0]) - - -def vthx(x): - return 0.01 - - -def vthy(x): - return 0.01 - - -def vthz(x): - return 0.01 - - -vvv = { - "vbulkx": vx, - "vbulky": vy, - "vbulkz": vz, - "vthx": vthx, - "vthy": vthy, - "vthz": vthz, -} - -ph.MaxwellianFluidModel( - bx=bx, - by=by, - bz=bz, - protons={"charge": 1, "density": density, **vvv, "init": {"seed": 1337}}, -) - -ph.ElectronModel(closure="isothermal", Te=0.12) - - -sim = ph.global_vars.sim - -timestamps = np.arange(0, sim.final_time + sim.time_step, 100 * sim.time_step) - - -for quantity in ["E", "B"]: - ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - - -for quantity in ["charge_density", "bulkVelocity"]: - ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) - -pops = [ - "protons", -] diff --git a/src/phare/phare_init_small.py b/src/phare/phare_init_small.py deleted file mode 100644 index f3f8e76e1..000000000 --- a/src/phare/phare_init_small.py +++ /dev/null @@ -1,112 +0,0 @@ -#!/usr/bin/env python3 - -import numpy as np - -import pyphare.pharein as ph - -# configure the simulation - -ph.Simulation( - smallest_patch_size=20, - largest_patch_size=20, - time_step_nbr=300, # number of time steps (not specified if time_step and final_time provided) - time_step=0.001, # simulation final time (not specified if time_step and time_step_nbr provided) - boundary_types="periodic", # boundary condition, string or tuple, length == len(cell) == len(dl) - cells=40, # integer or tuple length == dimension - dl=0.3, # mesh size of the root level, float or tuple - # max_nbr_levels=2, # (default=1) max nbr of levels in the AMR hierarchy - # refinement = "tagging", - refinement_boxes={"L0": {"B0": [(10,), (20,)]}}, - diag_options={ - "format": "phareh5", - "options": {"dir": "phare_outputs", "mode": "overwrite"}, - }, -) - - -def density(x): - return 1.0 - - -def by(x): - from pyphare.pharein.global_vars import sim - - L = sim.simulation_domain() - return 0.1 * np.cos(2 * np.pi * x / L[0]) - - -def bz(x): - from pyphare.pharein.global_vars import sim - - L = sim.simulation_domain() - return 0.1 * np.sin(2 * np.pi * x / L[0]) - - -def bx(x): - return 1.0 - - -def vx(x): - return 0.0 - - -def vy(x): - from pyphare.pharein.global_vars import sim - - L = sim.simulation_domain() - return 0.1 * np.cos(2 * np.pi * x / L[0]) - - -def vz(x): - from pyphare.pharein.global_vars import sim - - L = sim.simulation_domain() - return 0.1 * np.sin(2 * np.pi * x / L[0]) - - -def vthx(x): - return 0.01 - - -def vthy(x): - return 0.01 - - -def vthz(x): - return 0.01 - - -vvv = { - "vbulkx": vx, - "vbulky": vy, - "vbulkz": vz, - "vthx": vthx, - "vthy": vthy, - "vthz": vthz, -} - -ph.MaxwellianFluidModel( - bx=bx, - by=by, - bz=bz, - protons={"charge": 1, "density": density, **vvv, "init": {"seed": 1337}}, -) - -ph.ElectronModel(closure="isothermal", Te=0.12) - - -sim = ph.global_vars.sim - -timestamps = np.arange(0, sim.final_time + sim.time_step, 10 * sim.time_step) - - -for quantity in ["E", "B"]: - ph.ElectromagDiagnostics(quantity=quantity, write_timestamps=timestamps) - - -for quantity in ["charge_density", "bulkVelocity"]: - ph.FluidDiagnostics(quantity=quantity, write_timestamps=timestamps) - -pops = [ - "protons", -] diff --git a/tests/simulator/CMakeLists.txt b/tests/simulator/CMakeLists.txt index dde92f091..23da715ac 100644 --- a/tests/simulator/CMakeLists.txt +++ b/tests/simulator/CMakeLists.txt @@ -18,6 +18,7 @@ if(HighFive) phare_python3_exec(9 diagnostics test_diagnostics.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 phare_python3_exec(9 restarts test_restarts.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 phare_python3_exec(9 run test_run.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 + phare_python3_exec(9 vtk_diagnostics test_vtk_diagnostics.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 # phare_python3_exec(9 tagging test_tagging.py ${CMAKE_CURRENT_BINARY_DIR}) # serial or n = 2 if(testMPI) @@ -27,10 +28,11 @@ if(HighFive) # doesn't make sense in serial phare_mpi_python3_exec(9 3 load_balancing test_load_balancing.py ${CMAKE_CURRENT_BINARY_DIR}) - endif(testMPI) + phare_mpi_python3_exec(9 4 vtk_diagnostics test_vtk_diagnostics.py ${CMAKE_CURRENT_BINARY_DIR}) + endif(testMPI) - # phare_python3_exec(11 test_diagnostic_timestamps test_diagnostic_timestamps.py ${CMAKE_CURRENT_BINARY_DIR}) + phare_python3_exec(11 test_diagnostic_timestamps test_diagnostic_timestamps.py ${CMAKE_CURRENT_BINARY_DIR}) endif() diff --git a/tests/simulator/__init__.py b/tests/simulator/__init__.py index 3e131f81a..bb13dbf22 100644 --- a/tests/simulator/__init__.py +++ b/tests/simulator/__init__.py @@ -1,12 +1,12 @@ # # -import os import unittest import numpy as np from datetime import datetime import pyphare.pharein as ph +from pyphare.core.box import Box from pyphare.pharein import ElectronModel @@ -250,7 +250,7 @@ def unique_diag_dir_for_test_case(self, base_path, ndim, interp, post_path=""): from pyphare.cpp import cpp_lib cpp = cpp_lib() - return f"{base_path}/{self._testMethodName}/{cpp.mpi_size()}/{ndim}/{interp}/{post_path}" + return f"{base_path}/{self._testMethodName}/{cpp.mpi_size()}/{ndim}/{interp}{'/'+post_path if post_path else ''}" def clean_up_diags_dirs(self): from pyphare.cpp import cpp_lib diff --git a/tests/simulator/initialize/density_check.py b/tests/simulator/initialize/density_check.py index 7e474dedb..027ac6361 100644 --- a/tests/simulator/initialize/density_check.py +++ b/tests/simulator/initialize/density_check.py @@ -1,24 +1,16 @@ import os - -from pyphare.simulator.simulator import Simulator -from pyphare.pharesee.run import Run +import numpy as np +import matplotlib as mpl import pyphare.pharein as ph - +from pyphare.pharesee.run import Run from pyphare.pharein import global_vars -from tests.diagnostic import all_timestamps - +from pyphare.simulator.simulator import Simulator +from pyphare.pharesee.hierarchy import hierarchy_from from pyphare.pharein.diagnostics import FluidDiagnostics +from pyphare.pharesee.hierarchy import fromfunc -import matplotlib.pyplot as plt -import matplotlib as mpl -import numpy as np - -from pyphare.pharesee.hierarchy import hierarchy_from -from pyphare.pharesee.hierarchy.fromfunc import ions_mass_density_func1d -from pyphare.pharesee.hierarchy.fromfunc import ions_charge_density_func1d -from pyphare.pharesee.hierarchy.fromfunc import ions_mass_density_func2d -from pyphare.pharesee.hierarchy.fromfunc import ions_charge_density_func2d +from tests.diagnostic import all_timestamps mpl.use("Agg") @@ -219,6 +211,8 @@ def config_2d(): def main(): + import matplotlib.pyplot as plt + Simulator(config_1d()).run().reset() ph.global_vars.sim = None Simulator(config_2d()).run().reset() @@ -248,13 +242,13 @@ def assert_close_enough(h, H): H1 = hierarchy_from( hier=h1, - func=ions_mass_density_func1d, + func=fromfunc.ions_mass_density_func1d, masses=masses, densities=(densityMain_1d, densityBeam_1d), ) H2 = hierarchy_from( hier=h2, - func=ions_charge_density_func1d, + func=fromfunc.ions_charge_density_func1d, charges=charges, densities=(densityMain_1d, densityBeam_1d), ) @@ -283,13 +277,13 @@ def assert_close_enough(h, H): H1 = hierarchy_from( hier=h1, - func=ions_mass_density_func2d, + func=fromfunc.ions_mass_density_func2d, masses=masses, densities=(densityMain_2d, densityBeam_2d), ) H2 = hierarchy_from( hier=h2, - func=ions_charge_density_func2d, + func=fromfunc.ions_charge_density_func2d, charges=charges, densities=(densityMain_2d, densityBeam_2d), ) diff --git a/tests/simulator/refinement/test_2d_10_core.py b/tests/simulator/refinement/test_2d_10_core.py index 79c1aab5b..44ffb69d4 100644 --- a/tests/simulator/refinement/test_2d_10_core.py +++ b/tests/simulator/refinement/test_2d_10_core.py @@ -14,6 +14,14 @@ import pyphare.core.box as boxm from pyphare.simulator.simulator import Simulator, startMPI +from pyphare.cpp import cpp_lib +from tests.simulator.test_advance import AdvanceTestBase + +cpp = cpp_lib() +test = AdvanceTestBase(rethrow=True) # change to False for debugging images +L0_diags = "phare_outputs/test_x_homo_0" +L0L1_diags = "phare_outputs/test_x_homo_1" + ph.NO_GUI() @@ -120,16 +128,6 @@ def get_hier(path): return get_time(path) -from pyphare.cpp import cpp_lib - -from tests.simulator.test_advance import AdvanceTestBase - -cpp = cpp_lib() -test = AdvanceTestBase(rethrow=True) # change to False for debugging images -L0_diags = "phare_outputs/test_x_homo_0" -L0L1_diags = "phare_outputs/test_x_homo_1" - - def make_fig(hier, fig_name, ilvl, extra_collections=[]): if cpp.mpi_rank() == 0: l0_in_l1 = [boxm.refine(p.box, 2) for p in hier.level(0).patches] diff --git a/tests/simulator/test_restarts.py b/tests/simulator/test_restarts.py index 47eaa5c55..34e600a18 100644 --- a/tests/simulator/test_restarts.py +++ b/tests/simulator/test_restarts.py @@ -2,10 +2,10 @@ # import copy -import time import datetime import unittest import numpy as np +from time import sleep from pathlib import Path from datetime import timedelta @@ -349,7 +349,7 @@ def test_restarts_elapsed_time(self, ndim, interp, simInput, expected_num_levels # autodump false to ignore possible init dump simulator = Simulator(ph.global_vars.sim, auto_dump=False).initialize() - time.sleep(5) + sleep(5) simulator.advance().dump() # should trigger restart on "restart_idx" advance simulator.advance().dump() simulator.reset() @@ -431,7 +431,7 @@ def test_advanced_restarts_options(self): Dim / interp / etc are not relevant here """ ndim, interp = 1, 1 - print(f"test_advanced_restarts_options") + print("test_advanced_restarts_options") simput = copy.deepcopy( dup( @@ -456,10 +456,9 @@ def test_advanced_restarts_options(self): ph.global_vars.sim = None ph.global_vars.sim = ph.Simulation(**simput) - model = setup_model() + setup_model() Simulator(ph.global_vars.sim).run().reset() self.register_diag_dir_for_cleanup(local_out) - diag_dir0 = local_out simput["restart_options"]["restart_time"] = "auto" self.assertEqual(0.01, ph.restarts.restart_time(simput["restart_options"])) diff --git a/tests/simulator/test_run.py b/tests/simulator/test_run.py index 7d54c6fc6..c5153e15d 100644 --- a/tests/simulator/test_run.py +++ b/tests/simulator/test_run.py @@ -19,18 +19,15 @@ final_time = 0.05 time_step_nbr = int(final_time / time_step) 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(): +def config(diag_dir, diag_format): L = 0.5 sim = ph.Simulation( time_step=time_step, final_time=final_time, - cells=(40, 40), + cells=(20, 20), dl=(0.40, 0.40), refinement="tagging", max_nbr_levels=3, @@ -40,7 +37,7 @@ def config(): hyper_resistivity=0.002, resistivity=0.001, diag_options={ - "format": "phareh5", + "format": diag_format, "options": {"dir": diag_dir, "mode": "overwrite"}, }, ) @@ -163,60 +160,63 @@ def vthz(x, y): return sim -def plot_file_for_qty(qty, time): +def plot_file_for_qty(plot_dir, qty, time): return f"{plot_dir}/harris_{qty}_t{time}.png" def plot(diag_dir): run = Run(diag_dir) + plot_dir = Path(f"{diag_dir}_plots") + plot_dir.mkdir(parents=True, exist_ok=True) + pop_name = "protons" for time in timestamps: run.GetDivB(time).plot( - filename=plot_file_for_qty("divb", time), + filename=plot_file_for_qty(plot_dir, "divb", time), plot_patches=True, vmin=1e-11, vmax=2e-10, ) run.GetRanks(time).plot( - filename=plot_file_for_qty("Ranks", time), plot_patches=True + filename=plot_file_for_qty(plot_dir, "Ranks", time), plot_patches=True ) run.GetN(time, pop_name=pop_name).plot( - filename=plot_file_for_qty("N", time), plot_patches=True + filename=plot_file_for_qty(plot_dir, "N", time), plot_patches=True ) for c in ["x", "y", "z"]: run.GetB(time, all_primal=False).plot( - filename=plot_file_for_qty(f"b{c}", time), + filename=plot_file_for_qty(plot_dir, 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), + filename=plot_file_for_qty(plot_dir, "jz", time), qty="z", plot_patches=True, vmin=-2, vmax=2, ) run.GetPressure(time, pop_name=pop_name).plot( - filename=plot_file_for_qty(f"{pop_name}_Pxx", time), + filename=plot_file_for_qty(plot_dir, f"{pop_name}_Pxx", time), qty=pop_name + "_Pxx", plot_patches=True, ) run.GetPressure(time, pop_name=pop_name).plot( - filename=plot_file_for_qty(f"{pop_name}_Pzz", time), + filename=plot_file_for_qty(plot_dir, f"{pop_name}_Pzz", time), qty=pop_name + "_Pzz", plot_patches=True, ) run.GetPi(time).plot( - filename=plot_file_for_qty("Pxx", time), + filename=plot_file_for_qty(plot_dir, "Pxx", time), qty="Pxx", plot_patches=True, ) run.GetPi(time).plot( - filename=plot_file_for_qty("Pzz", time), + filename=plot_file_for_qty(plot_dir, "Pzz", time), qty="Pzz", plot_patches=True, ) + return plot_dir def assert_file_exists_with_size_at_least(file, size=10000): @@ -242,32 +242,44 @@ def tearDown(self): self.simulator = None ph.global_vars.sim = None - def test_run(self): - sim = config() + def _test_any_format(self, sim, diag_dir): self.register_diag_dir_for_cleanup(diag_dir) Simulator(sim).run().reset() - run = Run(diag_dir) + self.assertTrue(all(run.times("EM_B") == timestamps)) B = run.GetB(timestamps[-1], all_primal=False) self.assertTrue(B.levels()[0].patches[0].attrs) B = run.GetB(timestamps[-1]) self.assertTrue(B.levels()[0].patches[0].attrs) + def test_run_phareh5(self): + diag_dir = "phare_outputs/test_run_phareh5" + sim = config(diag_dir, "phareh5") + self._test_any_format(sim, diag_dir) + + # move to _test_any_format when vtkhdf supports divb etc if cpp.mpi_rank() == 0: - plot(diag_dir) + plot_dir = plot(diag_dir) for time in timestamps: for q in ["divb", "Ranks", "N", "jz"]: - assert_file_exists_with_size_at_least(plot_file_for_qty(q, time)) + assert_file_exists_with_size_at_least( + plot_file_for_qty(plot_dir, q, time) + ) for c in ["x", "y", "z"]: assert_file_exists_with_size_at_least( - plot_file_for_qty(f"b{c}", time) + plot_file_for_qty(plot_dir, f"b{c}", time) ) cpp.mpi_barrier() + def test_run_pharevtkhdf(self): + diag_dir = "phare_outputs/test_run_pharevtkhdf" + sim = config(diag_dir, "pharevtkhdf") + self._test_any_format(sim, diag_dir) + if __name__ == "__main__": import unittest diff --git a/tests/simulator/test_vtk_diagnostics.py b/tests/simulator/test_vtk_diagnostics.py new file mode 100644 index 000000000..00f18e4f3 --- /dev/null +++ b/tests/simulator/test_vtk_diagnostics.py @@ -0,0 +1,264 @@ +# + +import unittest +import itertools +import numpy as np +from copy import deepcopy +from ddt import data, ddt, unpack + +import pyphare.pharein as ph +from pyphare.cpp import cpp_lib +from pyphare.pharesee.run import Run + +from pyphare.core import phare_utilities as phut +from pyphare.simulator.simulator import startMPI +from pyphare.simulator.simulator import Simulator +from pyphare.pharesee.hierarchy import hierarchy_utils as hootils + +from tests.simulator import SimulatorTest +from tests.diagnostic import dump_all_diags + +cpp = cpp_lib() + +ndims = [2] # 3d todo / 1d not supported +interp_orders = [1] # , 2, 3 + + +def setup_model(sim, ppc=100): + L = 0.5 + + def density(x, y): + Ly = sim.simulation_domain()[1] + return ( + 0.4 + + 1.0 / np.cosh((y - Ly * 0.3) / L) ** 2 + + 1.0 / np.cosh((y - Ly * 0.7) / L) ** 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] + sigma = 1.0 + dB = 0.1 + + x0 = x - 0.5 * Lx + y1 = y - 0.3 * Ly + y2 = y - 0.7 * Ly + + dBy1 = 2 * dB * x0 * np.exp(-(x0**2 + y1**2) / (sigma) ** 2) + dBy2 = -2 * dB * x0 * np.exp(-(x0**2 + y2**2) / (sigma) ** 2) + + return dBy1 + dBy2 + + def bx(x, y): + Lx = sim.simulation_domain()[0] + Ly = sim.simulation_domain()[1] + sigma = 1.0 + dB = 0.1 + + x0 = x - 0.5 * Lx + y1 = y - 0.3 * Ly + y2 = y - 0.7 * Ly + + dBx1 = -2 * dB * y1 * np.exp(-(x0**2 + y1**2) / (sigma) ** 2) + dBx2 = 2 * dB * y2 * np.exp(-(x0**2 + y2**2) / (sigma) ** 2) + + v1 = -1 + v2 = 1.0 + return v1 + (v2 - v1) * (S(y, Ly * 0.3, L) - S(y, Ly * 0.7, L)) + dBx1 + dBx2 + + 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 = 0.7 + 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, + } + + model = ph.MaxwellianFluidModel( + bx=bx, + by=by, + bz=bz, + protons={ + "mass": 1, + "charge": 1, + "density": density, + **vvv, + "nbr_part_per_cell": ppc, + "init": {"seed": 1337}, + }, + alpha={ + "mass": 4, + "charge": 1, + "density": density, + **vvv, + "nbr_part_per_cell": ppc, + "init": {"seed": 2334}, + }, + ) + ph.ElectronModel(closure="isothermal", Te=0.12) + return model + + +out = "phare_outputs/vtk_diagnostic_test" +simArgs = { + "time_step_nbr": 1, + "final_time": 0.001, + "boundary_types": "periodic", + "cells": 40, + "dl": 0.3, + "diag_options": { + "format": "pharevtkhdf", + "options": {"dir": out, "mode": "overwrite"}, + }, +} + + +def permute(dic): + dic.update(simArgs.copy()) + return [ + dict( + ndim=ndim, + interp=interp_order, + simInput=deepcopy(dic), + ) + for ndim, interp_order in itertools.product(ndims, interp_orders) + ] + + +@ddt +class VTKDiagnosticsTest(SimulatorTest): + def __init__(self, *args, **kwargs): + super(VTKDiagnosticsTest, self).__init__(*args, **kwargs) + self.simulator = None + + def _run(self, ndim, interp, simInput, diag_dir="", **kwargs): + for key in ["cells", "dl", "boundary_types"]: + simInput[key] = list(phut.np_array_ify(simInput[key], ndim)) + local_out = self.unique_diag_dir_for_test_case( + f"{out}{'/'+diag_dir if diag_dir else ''}", ndim, interp + ) + self.register_diag_dir_for_cleanup(local_out) + simInput["diag_options"]["options"]["dir"] = local_out + simulation = ph.Simulation(**simInput) + self.assertTrue(len(simulation.cells) == ndim) + dump_all_diags(setup_model(simulation).populations) + Simulator(simulation).run().reset() + ph.global_vars.sim = None + return local_out + + @data(*permute({})) + @unpack + def test_dump_diags(self, ndim, interp, simInput): + print("test_dump_diags dim/interp:{}/{}".format(ndim, interp)) + + b0 = [[10 for i in range(ndim)], [19 for i in range(ndim)]] + simInput["refinement_boxes"] = {"L0": {"B0": b0}} + local_out = self._run(ndim, interp, simInput) + + try: + from pyphare.pharesee.phare_vtk import plot as plot_vtk + + plot_vtk(local_out + "/EM_B.vtkhdf", f"B{ndim}d.vtk.png") + plot_vtk(local_out + "/EM_E.vtkhdf", f"E{ndim}d.vtk.png") + except ModuleNotFoundError: + print("WARNING: vtk python module not found - cannot make plots") + + @data(*permute({})) + @unpack + def test_compare_l0_primal_to_phareh5(self, ndim, interp, simInput): + print("test_compare_l0_primal_to_phareh5 dim/interp:{}/{}".format(ndim, interp)) + + b0 = [[10 for i in range(ndim)], [19 for i in range(ndim)]] + simInput["refinement_boxes"] = {"L0": {"B0": b0}} + vtk_diags = self._run(ndim, interp, simInput, "test_vtk") + + simInput["diag_options"]["format"] = "phareh5" + phareh5_diags = self._run(ndim, interp, simInput, "test_h5") + + # not binary == with more than 2 cores + atol = 0 if cpp.mpi_size() <= 2 else 1e-16 + + for time in [0, simInput["final_time"]]: + # choose all primal value generally + phareh5_hier = Run(phareh5_diags).GetVi(time) + vtk_hier = Run(vtk_diags).GetVi(time) + + eqr = hootils.hierarchy_compare(phareh5_hier, vtk_hier, atol) + if not eqr: + print(eqr) + self.assertTrue(eqr) + + @data(*permute({})) + @unpack + def test_missing_level_case(self, ndim, interp, simInput): + simInput.update( + dict( + refinement="tagging", + max_nbr_levels=2, + tagging_threshold=0.99, # prevent level, + ) + ) + simInput["restart_options"] = dict( + dir="phare_outputs/vtk_padding_test", mode="overwrite", timestamps=[0.001] + ) + vtk_diags = self._run(ndim, interp, simInput, "test_vtk") + + simInput.update( + dict( + final_time=0.002, + tagging_threshold=0.01, # prefer level + ) + ) + + simInput["restart_options"]["timestamps"] = [] + simInput["restart_options"]["restart_time"] = 0.001 + del simInput["diag_options"]["options"]["mode"] # do not truncate diags + vtk_diags = self._run(ndim, interp, simInput, "test_vtk") + + hier0 = Run(vtk_diags).GetVi(0) + hier1 = Run(vtk_diags).GetVi(0.001) + hier2 = Run(vtk_diags).GetVi(0.002) + + self.assertTrue(1 not in hier0.levels()) + self.assertTrue(1 not in hier1.levels()) + self.assertTrue(1 in hier2.levels()) + + +if __name__ == "__main__": + startMPI() + unittest.main()