diff --git a/pyphare/pyphare/core/operators.py b/pyphare/pyphare/core/operators.py new file mode 100644 index 000000000..fcc20e5ad --- /dev/null +++ b/pyphare/pyphare/core/operators.py @@ -0,0 +1,147 @@ +import numpy as np + +from pyphare.pharesee.hierarchy import ScalarField, VectorField # TensorField +from pyphare.pharesee.hierarchy import compute_hier_from +from pyphare.pharesee.hierarchy import rename + + +def _compute_dot_product(patch_datas, **kwargs): + ref_name = next(iter(patch_datas.keys())) + + dset = ( + patch_datas["left_x"].dataset[:] * patch_datas["right_x"].dataset[:] + + patch_datas["left_y"].dataset[:] * patch_datas["right_y"].dataset[:] + + patch_datas["left_z"].dataset[:] * patch_datas["right_z"].dataset[:] + ) + + return ( + {"name": "value", "data": dset, "centering": patch_datas[ref_name].centerings}, + ) + + +def _compute_sqrt(patch_datas, **kwargs): + ref_name = next(iter(patch_datas.keys())) + + dset = np.sqrt(patch_datas["value"].dataset[:]) + + return ( + {"name": "value", "data": dset, "centering": patch_datas[ref_name].centerings}, + ) + + +def _compute_cross_product(patch_datas, **kwargs): + ref_name = next(iter(patch_datas.keys())) + + dset_x = ( + patch_datas["left_y"].dataset[:] * patch_datas["right_z"].dataset[:] + - patch_datas["left_z"].dataset[:] * patch_datas["right_y"].dataset[:] + ) + dset_y = ( + patch_datas["left_z"].dataset[:] * patch_datas["right_x"].dataset[:] + - patch_datas["left_x"].dataset[:] * patch_datas["right_z"].dataset[:] + ) + dset_z = ( + patch_datas["left_x"].dataset[:] * patch_datas["right_y"].dataset[:] + - patch_datas["left_y"].dataset[:] * patch_datas["right_x"].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_grad(patch_data, **kwargs): + ndim = patch_data["value"].box.ndim + nb_ghosts = kwargs["nb_ghosts"] + ds = patch_data["value"].dataset + + ds_shape = list(ds.shape) + + ds_x = np.full(ds_shape, np.nan) + ds_y = np.full(ds_shape, np.nan) + ds_z = np.full(ds_shape, np.nan) + + grad_ds = np.gradient(ds) + + if ndim == 2: + ds_x[nb_ghosts:-nb_ghosts, nb_ghosts:-nb_ghosts] = np.asarray( + grad_ds[0][nb_ghosts:-nb_ghosts, nb_ghosts:-nb_ghosts] + ) + ds_y[nb_ghosts:-nb_ghosts, nb_ghosts:-nb_ghosts] = np.asarray( + grad_ds[1][nb_ghosts:-nb_ghosts, nb_ghosts:-nb_ghosts] + ) + ds_z[nb_ghosts:-nb_ghosts, nb_ghosts:-nb_ghosts].fill( + 0.0 + ) # TODO at 2D, gradient is null in z dir + + else: + raise RuntimeError("dimension not yet implemented") + + return ( + {"name": "x", "data": ds_x, "centering": patch_data["value"].centerings}, + {"name": "y", "data": ds_y, "centering": patch_data["value"].centerings}, + {"name": "z", "data": ds_z, "centering": patch_data["value"].centerings}, + ) + + +def dot(hier_left, hier_right, **kwargs): + if isinstance(hier_left, VectorField) and isinstance(hier_right, VectorField): + names_left = ["left_x", "left_y", "left_z"] + names_right = ["right_x", "right_y", "right_z"] + + else: + raise RuntimeError("type of hierarchy not yet considered") + + hl = rename(hier_left, names_left) + hr = rename(hier_right, names_right) + + h = compute_hier_from( + _compute_dot_product, + (hl, hr), + ) + + return ScalarField(h) + + +def cross(hier_left, hier_right, **kwargs): + if isinstance(hier_left, VectorField) and isinstance(hier_right, VectorField): + names_left = ["left_x", "left_y", "left_z"] + names_right = ["right_x", "right_y", "right_z"] + + else: + raise RuntimeError("type of hierarchy not yet considered") + + hl = rename(hier_left, names_left) + hr = rename(hier_right, names_right) + + h = compute_hier_from( + _compute_cross_product, + (hl, hr), + ) + + return VectorField(h) + + +def sqrt(hier, **kwargs): + h = compute_hier_from( + _compute_sqrt, + hier, + ) + + return ScalarField(h) + + +def modulus(hier): + assert isinstance(hier, VectorField) + + return sqrt(dot(hier, hier)) + + +def grad(hier, **kwargs): + assert isinstance(hier, ScalarField) + nb_ghosts = list(hier.level(0).patches[0].patch_datas.values())[0].ghosts_nbr[0] + h = compute_hier_from(_compute_grad, hier, nb_ghosts=nb_ghosts) + + return VectorField(h) diff --git a/pyphare/pyphare/pharesee/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy.py index 2a4d1eb23..3c51da34b 100644 --- a/pyphare/pyphare/pharesee/hierarchy.py +++ b/pyphare/pyphare/pharesee/hierarchy.py @@ -563,6 +563,29 @@ def finest_part_data(hierarchy, time=None): class PatchHierarchy(object): """is a collection of patch levels""" + def flat_name(self, qty): + name_dict = { + "rho": "scalar", + "tags": "scalar", + "Bx": "x", + "Ex": "x", + "Fx": "x", + "Vx": "x", + "By": "y", + "Ey": "y", + "Fy": "y", + "Vy": "y", + "Bz": "z", + "Ez": "z", + "Fz": "z", + "Vz": "z", + } + + if qty in name_dict.keys(): + return name_dict[qty] + else: + raise ValueError("{} is not a valid quantity".format(qty)) + def __init__( self, patch_levels, @@ -572,10 +595,20 @@ def __init__( data_files=None, **kwargs, ): + times = time + if not isinstance(times, (list, tuple)): + times = listify(time) + + if not isinstance(patch_levels, (list, tuple)): + patch_levels = listify(patch_levels) + + assert len(times) == len(patch_levels) + self.patch_levels = patch_levels self.ndim = len(domain_box.lower) self.time_hier = {} - self.time_hier.update({self.format_timestamp(time): patch_levels}) + # self.time_hier.update({self.format_timestamp(time): patch_levels}) + self.time_hier.update({self.format_timestamp(t): pl for t, pl in zip(times, patch_levels)}) self.domain_box = domain_box self.refinement_ratio = refinement_ratio @@ -753,10 +786,10 @@ def global_min(self, qty, **kwargs): for patch in lvl.patches: pd = patch.patch_datas[qty] if first: - m = pd.dataset[:].min() + m = np.nanmin(pd.dataset[:]) first = False else: - m = min(m, pd.dataset[:].min()) + m = min(m, np.nanmin(pd.dataset[:])) return m @@ -767,10 +800,10 @@ def global_max(self, qty, **kwargs): for patch in lvl.patches: pd = patch.patch_datas[qty] if first: - m = pd.dataset[:].max() + m = np.nanmax(pd.dataset[:]) first = False else: - m = max(m, pd.dataset[:].max()) + m = max(m, np.nanmax(pd.dataset[:])) return m @@ -809,6 +842,9 @@ def __str__(self): s = s + "\n" return s + def get_names(self): + return list(self.levels()[0].patches[0].patch_datas.keys()) + def times(self): return np.sort(np.asarray(list(self.time_hier.keys()))) @@ -1088,6 +1124,394 @@ def dist_plot(self, **kwargs): return final, dp(final, **kwargs) + def __neg__(self): + names_self = self.get_names() + h = compute_hier_from(_compute_neg, self, names=names_self) + return VectorField(h) + + +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) + + +class ScalarField(PatchHierarchy): + + def __init__(self, hier): + renamed_hier = compute_hier_from(compute_rename, hier, new_names=("value",)) + patch_levels = renamed_hier.patch_levels + domain_box = renamed_hier.domain_box + refinement_ratio = renamed_hier.refinement_ratio + data_files = renamed_hier.data_files + + super().__init__(patch_levels, domain_box, refinement_ratio, renamed_hier.times(), data_files) + + def __add__(self, other): + assert isinstance(other, (ScalarField, int, float)) + h_self = rename(self, ["self_value"]) + + if isinstance(other, ScalarField): + h_other = rename(other, ["other_value"]) + h = compute_hier_from(self._compute_add, (h_self, h_other),) + elif isinstance(other, (int, float)): + h = compute_hier_from(self._compute_add, (h_self,), other=other) + pass + else: + raise RuntimeError("right operand not supported") + + return ScalarField(h) + + def __sub__(self, other): + assert isinstance(other, (ScalarField, int, float)) + h_self = rename(self, ["self_value"]) + + if isinstance(other, ScalarField): + h_other = rename(other, ["other_value"]) + h = compute_hier_from(self._compute_sub, (h_self, h_other),) + elif isinstance(other, (int, float)): + h = compute_hier_from(self._compute_sub, (h_self,), other=other) + else: + raise RuntimeError("right operand not supported") + + return ScalarField(h) + + def __mul__(self, other): + assert isinstance(other, (ScalarField, int, float)) + h_self = rename(self, ["self_value"]) + + if isinstance(other, ScalarField): + h_other = rename(other, ["other_value"]) + h = compute_hier_from(self._compute_mul, (h_self, h_other)) + elif isinstance(other, (int, float)): + h = compute_hier_from(self._compute_mul, (h_self,), other=other) + else: + raise RuntimeError("right operand not supported") + + return ScalarField(h) + + def __rmul__(self, other): + return self.__mul__(other) + + def __truediv__(self, other): + assert isinstance(other, (ScalarField, int, float)) + h_self = rename(self, ["self_value"]) + + if isinstance(other, ScalarField): + h_other = rename(other, ["other_value"]) + h = compute_hier_from(self._compute_truediv, (h_self, h_other)) + elif isinstance(other, (int, float)): + h = compute_hier_from(self._compute_truediv, (h_self,), other=other) + else: + raise RuntimeError("right operand not supported") + + return ScalarField(h) + + def __rtruediv__(self, other): + assert isinstance(other, (int, float)) + h_self = rename(self, ["self_value"]) + + h = compute_hier_from(self._compute_rtruediv, (h_self,), other=other) + + return ScalarField(h) + + def _compute_add(self, patch_datas, **kwargs): + ref_name = next(iter(patch_datas.keys())) + + if "other" in kwargs: + other = kwargs["other"] + dset_value = patch_datas["self_value"].dataset[:] + other + else: + dset_value = patch_datas["self_value"].dataset[:] +\ + patch_datas["other_value"].dataset[:] + + return ({"name": "value", + "data": dset_value, + "centering": patch_datas[ref_name].centerings},) + + def _compute_sub(self, patch_datas, **kwargs): + # subtract a scalar from the dataset of a scalarField + if "other" in kwargs: + other = kwargs["other"] + dset_value = patch_datas["self_value"].dataset[:] - other + # substraction of 2 scalarFields + else: + dset_value = patch_datas["self_value"].dataset[:] -\ + patch_datas["other_value"].dataset[:] + + return ({"name": "value", + "data": dset_value, + "centering": patch_datas["self_value"].centerings},) + + def _compute_mul(self, patch_datas, **kwargs): + # multiplication of a scalarField by a scalar + if "other" in kwargs: + other = kwargs["other"] + pd_attrs = [] + + for pd_name in patch_datas: + pd_attrs.append({"name": "value", + "data": other * patch_datas[pd_name].dataset[:], + "centering": patch_datas[pd_name].centerings, }) + # multiplication of 2 scalarField + else: + dset_value = patch_datas["self_value"].dataset[:] *\ + patch_datas["other_value"].dataset[:] + pd_attrs = ({"name": "value", + "data": dset_value, + "centering": patch_datas["self_value"].centerings},) + + return tuple(pd_attrs) + + def _compute_truediv(self, patch_datas, **kwargs): + # multiplication of a scalarField by a scalar + if "other" in kwargs: + other = kwargs["other"] + pd_attrs = [] + + for pd_name in patch_datas: + pd_attrs.append({"name": "value", + "data": patch_datas[pd_name].dataset[:] / other, + "centering": patch_datas[pd_name].centerings, }) + # multiplication of 2 scalarField + else: + dset_value = patch_datas["self_value"].dataset[:] /\ + patch_datas["other_value"].dataset[:] + pd_attrs = ({"name": "value", + "data": dset_value, + "centering": patch_datas["self_value"].centerings},) + + return tuple(pd_attrs) + + def _compute_rtruediv(self, patch_datas, **kwargs): + # Scalar divided by a scalarField + other = kwargs["other"] + pd_attrs = [] + + for pd_name in patch_datas: + pd_attrs.append({"name": "value", + "data": other / patch_datas[pd_name].dataset[:], + "centering": patch_datas[pd_name].centerings, }) + + return tuple(pd_attrs) + + +class VectorField(PatchHierarchy): + + def __init__(self, hier): + renamed_hier = compute_hier_from(compute_rename, hier, new_names=("x", "y", "z")) + patch_levels = renamed_hier.patch_levels + domain_box = renamed_hier.domain_box + refinement_ratio = renamed_hier.refinement_ratio + data_files = renamed_hier.data_files + + self.names = ["x", "y", "z"] + + super().__init__(patch_levels, domain_box, refinement_ratio, renamed_hier.times(), data_files) + + def __mul__(self, other): + assert isinstance(other, (int, float)) + h = compute_hier_from(_compute_mul, self, names=["x", "y", "z"], other=other) + return VectorField(h) + + def __rmul__(self, other): + return self.__mul__(other) + + def __add__(self, other): + names_self_kept = self.get_names() + names_other_kept = other.get_names() + + if isinstance(other, VectorField): + names_self = ["self_x", "self_y", "self_z"] + names_other = ["other_x", "other_y", "other_z"] + else: + raise RuntimeError("type of hierarchy not yet considered") + + h_self = rename(self, names_self) + h_other = rename(other, names_other) + + h = compute_hier_from(_compute_add, (h_self, h_other),) + + self = rename(h_self, names_self_kept) # needed ? + other = rename(h_other, names_other_kept) + + return VectorField(h) + + def __sub__(self, other): + names_self_kept = self.get_names() + names_other_kept = other.get_names() + + if isinstance(other, VectorField): + names_self = ["self_x", "self_y", "self_z"] + names_other = ["other_x", "other_y", "other_z"] + else: + raise RuntimeError("type of hierarchy not yet considered") + + h_self = rename(self, names_self) + h_other = rename(other, names_other) + + h = compute_hier_from(_compute_sub, (h_self, h_other),) + + self = rename(h_self, names_self_kept) + other = rename(h_other, names_other_kept) + + return VectorField(h) + + def __truediv__(self, other): + + if not isinstance(other, (ScalarField, int, float)): + raise RuntimeError("type of operand not considered") + + if isinstance(other, ScalarField): + return VectorField(compute_hier_from(_compute_truediv, + (self, other), + res_names=("x", "y", "z"))) + elif isinstance(other, (int, float)): + return VectorField(compute_hier_from(_compute_scalardiv, + (self,), + res_names=("x", "y", "z"), + scalar=other)) + + +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["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) + + +def _compute_rename(patch_datas, **kwargs): + new_names = kwargs["names"] + pd_attrs = [] + + for name, new_name in zip(patch_datas.keys(), new_names): + pd_attrs.append( + { + "name": new_name, + "data": patch_datas[name].dataset, + "centering": patch_datas[name].centerings, + } + ) + + return tuple(pd_attrs) + + +def rename(hierarchy, names): + return compute_hier_from(_compute_rename, hierarchy, names=names) + def amr_grid(hierarchy, time): """returns a non-uniform contiguous primal grid @@ -1274,11 +1698,11 @@ def are_compatible_hierarchies(hierarchies): return True -def extract_patchdatas(hierarchies, ilvl, ipatch): +def extract_patchdatas(hierarchies, ilvl, ipatch, t): """ returns a dict {patchdata_name:patchdata} from a list of hierarchies for patch ipatch at level ilvl """ - patches = [h.patch_levels[ilvl].patches[ipatch] for h in hierarchies] + patches = [h.level(ilvl,t).patches[ipatch] for h in hierarchies] patch_datas = { pdname: pd for p in patches for pdname, pd in list(p.patch_datas.items()) } @@ -1294,14 +1718,14 @@ def new_patchdatas_from(compute, patchdatas, layout, **kwargs): return new_patch_datas -def new_patches_from(compute, hierarchies, ilvl, **kwargs): +def new_patches_from(compute, hierarchies, ilvl,t, **kwargs): reference_hier = hierarchies[0] new_patches = [] - patch_nbr = len(reference_hier.patch_levels[ilvl].patches) + patch_nbr = len(reference_hier.level(ilvl, t).patches) for ip in range(patch_nbr): - current_patch = reference_hier.patch_levels[ilvl].patches[ip] + current_patch = reference_hier.level(ilvl, t).patches[ip] layout = current_patch.layout - patch_datas = extract_patchdatas(hierarchies, ilvl, ip) + patch_datas = extract_patchdatas(hierarchies, ilvl, ip, t) new_patch_datas = new_patchdatas_from(compute, patch_datas, layout, **kwargs) new_patches.append(Patch(new_patch_datas, current_patch.id)) return new_patches @@ -1317,11 +1741,12 @@ def compute_hier_from(compute, hierarchies, **kwargs): hierarchies = listify(hierarchies) reference_hier = hierarchies[0] domain_box = reference_hier.domain_box - patch_levels = {} - for ilvl in range(reference_hier.levelNbr()): - patch_levels[ilvl] = PatchLevel( - ilvl, new_patches_from(compute, hierarchies, ilvl, **kwargs) - ) + patch_levels = [{}]*len(reference_hier.times()) + for it, t in enumerate(reference_hier.times()): + for ilvl in range(reference_hier.levelNbr()): + patch_levels[it][ilvl] = PatchLevel( + ilvl, new_patches_from(compute, hierarchies, ilvl, t, **kwargs) + ) assert len(reference_hier.time_hier) == 1 # only single time hierarchies now t = list(reference_hier.time_hier.keys())[0] diff --git a/pyphare/pyphare/pharesee/run.py b/pyphare/pyphare/pharesee/run.py index 988479df9..96fbdf429 100644 --- a/pyphare/pyphare/pharesee/run.py +++ b/pyphare/pyphare/pharesee/run.py @@ -1,10 +1,15 @@ import os - import numpy as np -from pyphare.pharesee.hierarchy import compute_hier_from -from pyphare.pharesee.hierarchy import compute_hier_from_ -from .hierarchy import flat_finest_field, hierarchy_from +from pyphare.pharesee.hierarchy import ( + compute_hier_from, + flat_finest_field, + hierarchy_from, + ScalarField, + VectorField, +) + +from pyphare.core.gridlayout import yee_centering def _current1d(by, bz, xby, xbz): @@ -130,6 +135,234 @@ def _compute_divB(patchdatas, **kwargs): raise RuntimeError("dimension not implemented") +def _ppp_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (primal,primal,primal) to (primal,primal,primal) + centering that is the centering of moments on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, _, _ = _inner_slices(nb_ghosts) + + inner_all = tuple([inner] * ndim) + return inner_all, (inner_all,) + + +def _pdd_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (dual,primal,primal) to (primal,primal,primal) + centering that is the centering of Bx on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, inner_shift_left, inner_shift_right = _inner_slices(nb_ghosts) + + if ndim == 1: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_all,) + elif ndim == 2: + inner_all = tuple([inner] * ndim) + return inner_all, ((inner, inner_shift_left), + (inner, inner_shift_right)) + else: + raise RuntimeError("dimension not yet implemented") + + +def _dpd_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (dual,primal,primal) to (primal,primal,primal) + centering that is the centering of By on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, inner_shift_left, inner_shift_right = _inner_slices(nb_ghosts) + + if ndim == 1: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_shift_left, + inner_shift_right) + elif ndim == 2: + inner_all = tuple([inner] * ndim) + return inner_all, ((inner_shift_left, inner), + (inner_shift_right, inner)) + else: + raise RuntimeError("dimension not yet implemented") + + +def _ddp_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (dual,primal,primal) to (primal,primal,primal) + centering that is the centering of Bz on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, inner_shift_left, inner_shift_right = _inner_slices(nb_ghosts) + + if ndim == 1: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_shift_left, + inner_shift_right) + elif ndim == 2: + inner_all = tuple([inner] * ndim) + return inner_all, ((inner_shift_left, inner_shift_left), + (inner_shift_left, inner_shift_right), + (inner_shift_right, inner_shift_left), + (inner_shift_right, inner_shift_right)) + else: + raise RuntimeError("dimension not yet implemented") + + +def _dpp_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (dual,primal,primal) to (primal,primal,primal) + centering that is the centering of Ex on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, inner_shift_left, inner_shift_right = _inner_slices(nb_ghosts) + + if ndim == 1: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_shift_left, + inner_shift_right) + elif ndim == 2: + inner_all = tuple([inner] * ndim) + return inner_all, ((inner_shift_left, inner), + (inner_shift_right, inner)) + else: + raise RuntimeError("dimension not yet implemented") + + +def _pdp_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (dual,primal,primal) to (primal,primal,primal) + centering that is the centering of Ey on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, inner_shift_left, inner_shift_right = _inner_slices(nb_ghosts) + + if ndim == 1: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_all,) + elif ndim == 2: + inner_all = tuple([inner] * ndim) + return inner_all, ((inner, inner_shift_left), + (inner, inner_shift_right)) + else: + raise RuntimeError("dimension not yet implemented") + + +def _ppd_to_ppp_domain_slicing(**kwargs): + """ + return the slicing for (dual,primal,primal) to (primal,primal,primal) + centering that is the centering of Ez on a Yee grid + """ + + nb_ghosts = kwargs["nb_ghosts"] + ndim = kwargs["ndim"] + + inner, _, _ = _inner_slices(nb_ghosts) + + if ndim == 1: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_all,) + elif ndim == 2: + inner_all = tuple([inner] * ndim) + return inner_all, (inner_all,) + else: + raise RuntimeError("dimension not yet implemented") + + +slices_to_primal_ = {"primal_primal_primal": _ppp_to_ppp_domain_slicing, + "primal_dual_dual": _pdd_to_ppp_domain_slicing, + "dual_primal_dual": _dpd_to_ppp_domain_slicing, + "dual_dual_primal": _ddp_to_ppp_domain_slicing, + "dual_primal_primal": _dpp_to_ppp_domain_slicing, + "primal_dual_primal": _pdp_to_ppp_domain_slicing, + "primal_primal_dual": _ppd_to_ppp_domain_slicing, + } + + +def merge_centerings(pdname): + from pyphare.core.gridlayout import directions + return "_".join([yee_centering[d][pdname] for d in directions]) + + +def slices_to_primal(pdname, **kwargs): + return slices_to_primal_[merge_centerings(pdname)](**kwargs) + + +def _compute_to_primal(patchdatas, **kwargs): + """ + datasets have NaN in their ghosts... might need to be properly filled + with their neighbors already properly projected on primal + """ + + reference_name = next(iter(kwargs.values())) + reference_pd = patchdatas[reference_name] + nb_ghosts = reference_pd.layout.nbrGhosts( + reference_pd.layout.interp_order, "primal" + ) + ndim = reference_pd.box.ndim + + centerings = ["primal"] * ndim + + pd_attrs = [] + for name, pd_name in kwargs.items(): + pd = patchdatas[pd_name] + + ds = pd.dataset + + ds_shape = list(ds.shape) + for i in range(ndim): + if pd.centerings[i] == "dual": + ds_shape[i] += 1 + + # should be something else than nan values when the ghosts cells + # will be filled with correct values coming from the neighbors + ds_all_primal = np.full(ds_shape, np.nan) + ds_ = np.zeros(ds_shape) + + # inner is the slice containing the points that are updated + # in the all_primal dataset + # chunks is a tupls of all the slices coming from the initial dataset + # that are needed to calculate the average for the all_primal dataset + inner, chunks = slices_to_primal(pd_name, + nb_ghosts=nb_ghosts, + ndim=ndim) + + for chunk in chunks: + ds_[inner] = np.add(ds_[inner], ds[chunk] / len(chunks)) + ds_all_primal[inner] = ds_[inner] + + pd_attrs.append({"name": name, + "data": ds_all_primal, + "centering": centerings}) + + return tuple(pd_attrs) + + +def _inner_slices(nb_ghosts): + inner = slice(nb_ghosts, -nb_ghosts) + inner_shift_left = slice(nb_ghosts - 1, -nb_ghosts) + inner_shift_right = slice(nb_ghosts, -nb_ghosts + 1) + + return inner, inner_shift_left, inner_shift_right + + def _get_rank(patchdatas, **kwargs): """ make a field dataset cell centered coding the MPI rank @@ -137,7 +370,7 @@ def _get_rank(patchdatas, **kwargs): """ from pyphare.core.box import grow - reference_pd = patchdatas["Bx"] # take Bx as a reference, but could be any other + reference_pd = patchdatas["Bx"] # Bx as a ref, but could be any other ndim = reference_pd.box.ndim pid = kwargs["id"] @@ -321,35 +554,44 @@ def _get(self, hierarchy, time, merged, interp): def GetTags(self, time, merged=False): hier = self._get_hierarchy(time, "tags.h5") - return self._get(hier, time, merged, "nearest") + return ScalarField(self._get(hier, time, merged, "nearest")) + # return self._get(hier, time, merged, "nearest") - def GetB(self, time, merged=False, interp="nearest"): + def GetB(self, time, merged=False, interp="nearest", all_primal=True): hier = self._get_hierarchy(time, "EM_B.h5") - return self._get(hier, time, merged, interp) + if not all_primal: + return self._get(hier, time, merged, interp) + + h = compute_hier_from(_compute_to_primal, hier, x="Bx", y="By", z="Bz") + return VectorField(h) - def GetE(self, time, merged=False, interp="nearest"): + def GetE(self, time, merged=False, interp="nearest", all_primal=True): hier = self._get_hierarchy(time, "EM_E.h5") - return self._get(hier, time, merged, interp) + if not all_primal: + return self._get(hier, time, merged, interp) + + h = compute_hier_from(_compute_to_primal, hier, x="Ex", y="Ey", z="Ez") + return VectorField(h) def GetMassDensity(self, time, merged=False, interp="nearest"): hier = self._get_hierarchy(time, "ions_mass_density.h5") - return self._get(hier, time, merged, interp) + return ScalarField(self._get(hier, time, merged, interp)) def GetNi(self, time, merged=False, interp="nearest"): hier = self._get_hierarchy(time, "ions_density.h5") - return self._get(hier, time, merged, interp) + return ScalarField(self._get(hier, time, merged, interp)) def GetN(self, time, pop_name, merged=False, interp="nearest"): hier = self._get_hierarchy(time, f"ions_pop_{pop_name}_density.h5") - return self._get(hier, time, merged, interp) + return ScalarField(self._get(hier, time, merged, interp)) def GetVi(self, time, merged=False, interp="nearest"): hier = self._get_hierarchy(time, "ions_bulkVelocity.h5") - return self._get(hier, time, merged, interp) + return VectorField(self._get(hier, time, merged, interp)) def GetFlux(self, time, pop_name, merged=False, interp="nearest"): hier = self._get_hierarchy(time, f"ions_pop_{pop_name}_flux.h5") - return self._get(hier, time, merged, interp) + return VectorField(self._get(hier, time, merged, interp)) def GetPressure(self, time, pop_name, merged=False, interp="nearest"): M = self._get_hierarchy(time, f"ions_pop_{pop_name}_momentum_tensor.h5") @@ -361,24 +603,38 @@ def GetPressure(self, time, pop_name, merged=False, interp="nearest"): popname=pop_name, mass=self.GetMass(pop_name), ) - return self._get(P, time, merged, interp) + return self._get(P, time, merged, interp) # should later be a TensorField def GetPi(self, time, merged=False, interp="nearest"): M = self._get_hierarchy(time, f"ions_momentum_tensor.h5") massDensity = self.GetMassDensity(time) Vi = self._get_hierarchy(time, f"ions_bulkVelocity.h5") Pi = compute_hier_from(_compute_pressure, (M, massDensity, Vi)) - return self._get(Pi, time, merged, interp) + return self._get(Pi, time, merged, interp) # should later be a TensorField - def GetJ(self, time, merged=False, interp="nearest"): - B = self.GetB(time) + def GetPe(self, time, merged=False, interp="nearest", all_primal=True): + hier = self._get_hierarchy(time, "ions_density.h5") + + Te = hier.sim.electrons.closure.Te + + if not all_primal: + return Te * self._get(hier, time, merged, interp) + + h = compute_hier_from(_compute_to_primal, hier, scalar="rho") + return ScalarField(h)*Te + + def GetJ(self, time, merged=False, interp="nearest", all_primal=True): + B = self.GetB(time, all_primal=False) J = compute_hier_from(_compute_current, B) - return self._get(J, time, merged, interp) + if not all_primal: + return self._get(J, time, merged, interp) + h = compute_hier_from(_compute_to_primal, J, x="Jx", y="Jy", z="Jz") + return VectorField(h) def GetDivB(self, time, merged=False, interp="nearest"): - B = self.GetB(time) + B = self.GetB(time, all_primal=False) db = compute_hier_from(_compute_divB, B) - return self._get(db, time, merged, interp) + return ScalarField(self._get(db, time, merged, interp)) def GetRanks(self, time, merged=False, interp="nearest"): """ @@ -389,7 +645,7 @@ def GetRanks(self, time, merged=False, interp="nearest"): """ B = self.GetB(time) ranks = compute_hier_from(_get_rank, B) - return self._get(ranks, time, merged, interp) + return ScalarField(self._get(ranks, time, merged, interp)) def GetParticles(self, time, pop_name, hier=None): def filename(name): @@ -458,17 +714,24 @@ def GetDl(self, level="finest", time=None): return root_cell_width / fac - def GetAllAvailableQties(self, time, pops): + def GetAllAvailableQties(self, time, pops, all_primal=True): + assert not all_primal assert self.single_hier_for_all_quantities == True # can't work otherwise - self.GetParticles(time, pops) - self.GetB(time) - self.GetE(time) - self.GetNi(time) - self.GetVi(time) + def _try(fn, *args, **kwargs): + try: + fn(*args, **kwargs) + except: + ... # file not found + + _try(self.GetParticles, time, pops) + _try(self.GetB, time, all_primal=all_primal) + _try(self.GetE, time, all_primal=all_primal) + _try(self.GetNi, time, all_primal=all_primal) + _try(self.GetVi, time) for pop in pops: - self.GetFlux(time, pop) - self.GetN(time, pop) + _try(self.GetFlux, time, pop) + _try(self.GetN, time, pop) return self.hier diff --git a/subprojects/cppdict b/subprojects/cppdict deleted file mode 160000 index 6b02756b9..000000000 --- a/subprojects/cppdict +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 6b02756b9646811c65b5df83497b82ba75f1270e diff --git a/tests/core/data/gridlayout/test_laplacian.py b/tests/core/data/gridlayout/test_laplacian.py index 91c689475..1ef82816d 100644 --- a/tests/core/data/gridlayout/test_laplacian.py +++ b/tests/core/data/gridlayout/test_laplacian.py @@ -210,9 +210,7 @@ def test_laplacian_yee2D(path): Jz = np.tensordot(np.sinh(0.2 * x_primal), np.cosh(0.2 * y_primal), axes=0) Jx_x[psi_d_X : pei_d_X + 1, :] = ( - Jx[ - psi_d_X + 1 : pei_d_X + 2 :, - ] + Jx[psi_d_X + 1 : pei_d_X + 2 :,] - 2.0 * Jx[psi_d_X : pei_d_X + 1, :] + Jx[psi_d_X - 1 : pei_d_X, :] ) / (tv.meshSize[0] * tv.meshSize[0]) diff --git a/tests/simulator/test_restarts.py b/tests/simulator/test_restarts.py index f847ac56c..faf482347 100644 --- a/tests/simulator/test_restarts.py +++ b/tests/simulator/test_restarts.py @@ -166,8 +166,8 @@ def count_levels_and_patches(qty): run0 = Run(diag_dir0, single_hier_for_all_quantities=True) run1 = Run(diag_dir1, single_hier_for_all_quantities=True) - datahier0 = run0.GetAllAvailableQties(time, pops) - datahier1 = run1.GetAllAvailableQties(time, pops) + datahier0 = run0.GetAllAvailableQties(time, pops, all_primal=False) + datahier1 = run1.GetAllAvailableQties(time, pops, all_primal=False) self.assertEqual( datahier0.level(0).patches[0].patch_datas.keys(), @@ -208,7 +208,9 @@ def count_levels_and_patches(qty): np.testing.assert_equal(pd0.dataset[:], pd1.dataset[:]) checks += 1 - n_levels, n_patches = count_levels_and_patches(run0.GetB(time)) + n_levels, n_patches = count_levels_and_patches( + run0.GetB(time, all_primal=False) + ) self.assertEqual(n_levels, expected_num_levels) self.assertGreaterEqual(n_patches, n_levels) # at least one patch per level self.assertEqual(checks, n_quantities_per_patch * n_patches)