can set a patchData w. an all_primal centering

pyphare/pyphare/core/operators.py

@@ -0,0 +1,150 @@
+import numpy as np
+# from pyphare.pharesee.hierarchy import PatchLevel, Patch, FieldData
+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": 'scalar', "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["scalar"].dataset[:])
+
+    return ({"name": 'scalar', "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["scalar"].box.ndim
+    nb_ghosts = kwargs["nb_ghosts"]
+    ds = patch_data["scalar"].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.)  # 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["scalar"].centerings},
+            {"name": 'y', "data": ds_y, "centering": patch_data["scalar"].centerings},
+            {"name": 'z', "data": ds_z, "centering": patch_data["scalar"].centerings},
+            )
+
+
+def dot(hier_left, hier_right, **kwargs):
+    names_left_kept = hier_left.get_names()
+    names_right_kept = hier_right.get_names()
+
+    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),)
+
+    hier_left = rename(hl, names_left_kept)
+    hier_right = rename(hr, names_right_kept)
+
+    return ScalarField(h.patch_levels, h.domain_box,
+                       refinement_ratio=h.refinement_ratio,
+                       time=h.times()[0],
+                       data_files=h.data_files)
+
+
+def cross(hier_left, hier_right, **kwargs):
+    names_left_kept = hier_left.get_names()
+    names_right_kept = hier_right.get_names()
+
+    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),)
+
+    hier_left = rename(hl, names_left_kept)
+    hier_right = rename(hr, names_right_kept)
+
+    return VectorField(h.patch_levels, h.domain_box,
+                       refinement_ratio=h.refinement_ratio,
+                       time=h.times()[0],
+                       data_files=h.data_files)
+
+
+def sqrt(hier, **kwargs):
+    h = compute_hier_from(_compute_sqrt, hier,)
+
+    return ScalarField(h.patch_levels, h.domain_box,
+                       refinement_ratio=h.refinement_ratio,
+                       time=h.times()[0],
+                       data_files=h.data_files)
+
+
+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)
+
+    # TODO the plot of a grad displays only 1 patch if vmin and vmax are not specified... any idea why ?
+    return VectorField(h.patch_levels, h.domain_box,
+                       refinement_ratio=h.refinement_ratio,
+                       time=h.times()[0],
+                       data_files=h.data_files)

pyphare/pyphare/pharesee/hierarchy.py

@@ -563,6 +563,18 @@ def finest_part_data(hierarchy, time=None):
 class PatchHierarchy(object):
     """is a collection of patch levels"""
 
+    def flat_name(self, qty):
+        if qty in ['rho', 'tags']:
+            return 'scalar'
+        elif qty in ['Bx', 'Ex', 'Fx', 'Vx']:
+            return 'x'
+        elif qty in ['By', 'Ey', 'Fy', 'Vy']:
+            return 'y'
+        elif qty in ['Bz', 'Ez', 'Fz', 'Vz']:
+            return 'z'
+        else:
+            raise ValueError("{} is not a valid quantity".format(qty))
+
     def __init__(
         self,
         patch_levels,
@@ -790,6 +802,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())))
 
@@ -1069,6 +1084,209 @@ 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.patch_levels, h.domain_box,
+                           refinement_ratio=h.refinement_ratio,
+                           time=h.times()[0],
+                           data_files=h.data_files)
+
+
+class ScalarField(PatchHierarchy):
+
+    name = 'scalar'
+
+    def __init__(self, patch_levels, domain_box, **kwargs):
+        refinement_ratio = kwargs.get("refinement_ratio", 2)
+        time = kwargs.get("time", .0)
+        data_files = kwargs.get("data_files", None)
+
+        super().__init__(patch_levels, domain_box, refinement_ratio, time, data_files)
+
+    def __mul__(self, other):
+        assert isinstance(other, (int, float))
+
+        h = compute_hier_from(_compute_mul, self, names=self.get_names(), other=other)
+
+        return ScalarField(h.patch_levels, h.domain_box,
+                           refinement_ratio=h.refinement_ratio,
+                           time=h.times()[0],
+                           data_files=h.data_files)
+
+    def __rmul__(self, other):
+        return self.__mul__(other)
+
+
+class VectorField(PatchHierarchy):
+
+    def __init__(self, patch_levels, domain_box, **kwargs):
+        refinement_ratio = kwargs.get("refinement_ratio", 2)
+        time = kwargs.get("time", .0)
+        data_files = kwargs.get("data_files", None)
+
+        self.names = ['x', 'y', 'z']
+
+        super().__init__(patch_levels, domain_box, refinement_ratio, time, data_files)
+
+    def __mul__(self, other):
+        assert isinstance(other, (int, float))
+
+        h = compute_hier_from(_compute_mul, self, names=self.get_names(), other=other)
+
+        return VectorField(h.patch_levels, h.domain_box,
+                           refinement_ratio=h.refinement_ratio,
+                           time=h.times()[0],
+                           data_files=h.data_files)
+
+    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)
+        other = rename(h_other, names_other_kept)
+
+        return VectorField(h.patch_levels, h.domain_box,
+                           refinement_ratio=h.refinement_ratio,
+                           time=h.times()[0],
+                           data_files=h.data_files)
+
+    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.patch_levels, h.domain_box,
+                           refinement_ratio=h.refinement_ratio,
+                           time=h.times()[0],
+                           data_files=h.data_files)
+
+    def __truediv__(self, other):
+
+        if isinstance(other, ScalarField):
+            names = self.get_names()
+        else:
+            raise RuntimeError("type of hierarchy not yet considered")
+
+        h = compute_hier_from(_compute_truediv, (self, other), names=names)
+
+        return VectorField(h.patch_levels, h.domain_box,
+                           refinement_ratio=h.refinement_ratio,
+                           time=h.times()[0],
+                           data_files=h.data_files)
+
+
+def _compute_mul(patch_datas, **kwargs):
+    names = kwargs["names"]
+    other = kwargs["other"]
+    pd_attrs = []
+
+    for name in names:
+        pd_attrs.append({"name": name,
+                         "data": other*patch_datas[name].dataset,
+                         "centering": patch_datas[name].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["names"]
+    pd_attrs = []
+
+    for name in names:
+        pd_attrs.append({"name": name,
+                         "data": patch_datas[name].dataset/patch_datas["scalar"].dataset,
+                         "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