ufuncing

pyphare/pyphare/core/box.py

@@ -207,7 +207,7 @@ def amr_to_local(box, ref_box):
 
 
 def select(data, box):
-    return data[tuple([slice(l, u + 1) for l, u in zip(box.lower, box.upper)])]
+    return data[tuple([slice(lo, up + 1) for lo, up in zip(box.lower, box.upper)])]
 
 
 class DataSelector:

pyphare/pyphare/core/gridlayout.py

@@ -44,6 +44,10 @@
         "Pyz": "primal",
         "Pzz": "primal",
         "tags": "dual",
+        "value": "primal",
+        "x": "primal",
+        "y": "primal",
+        "z": "primal",
     },
     "y": {
         "Bx": "dual",
@@ -79,6 +83,10 @@
         "Pyz": "primal",
         "Pzz": "primal",
         "tags": "dual",
+        "value": "primal",
+        "x": "primal",
+        "y": "primal",
+        "z": "primal",
     },
     "z": {
         "Bx": "dual",
@@ -114,6 +122,10 @@
         "Pyz": "primal",
         "Pzz": "primal",
         "tags": "dual",
+        "value": "primal",
+        "x": "primal",
+        "y": "primal",
+        "z": "primal",
     },
 }
 yee_centering_lower = {

pyphare/pyphare/core/operators.py

@@ -96,7 +96,7 @@ def dot(hier_left, hier_right, **kwargs):
         (hl, hr),
     )
 
-    return ScalarField(h)
+    return ScalarField.FROM(h)
 
 
 def cross(hier_left, hier_right, **kwargs):
@@ -115,7 +115,7 @@ def cross(hier_left, hier_right, **kwargs):
         (hl, hr),
     )
 
-    return VectorField(h)
+    return VectorField.FROM(h)
 
 
 def sqrt(hier, **kwargs):
@@ -124,7 +124,7 @@ def sqrt(hier, **kwargs):
         hier,
     )
 
-    return ScalarField(h)
+    return ScalarField.FROM(h)
 
 
 def modulus(hier):
@@ -138,4 +138,4 @@ def grad(hier, **kwargs):
     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)
+    return VectorField.FROM(h)

pyphare/pyphare/pharesee/hierarchy/hierarchy.py

@@ -16,6 +16,40 @@ def format_timestamp(timestamp):
     return "{:.10f}".format(timestamp)
 
 
+class IndexHierarchy:
+    def __init__(self, time, hier, indexes):
+        self.time = time
+        self.hier = hier
+        self.indexes = indexes
+
+    def __str__(self):
+        s = "IndexHierarchy: \n"
+        s += "Time {}\n".format(self.time)
+        for ilvl, lvl in self.indexes.items():
+            s += "Level {}\n".format(ilvl)
+            for ip, qty_indexes in enumerate(lvl):
+                for qty_name, indexes in qty_indexes.items():
+                    s += f"    P{ip} {type} {qty_name} {indexes} \n"
+        return s
+
+
+class ValueHierarchy:
+    def __init__(self, time, hier, values):
+        self.time = time
+        self.hier = hier
+        self.values = values
+
+    def __str__(self):
+        s = "ValueHierarchy: \n"
+        s += "Time {}\n".format(self.time)
+        for ilvl, lvl in self.values.items():
+            s += "Level {}\n".format(ilvl)
+            for ip, qty_values in enumerate(lvl):
+                for qty_name, values in qty_values.items():
+                    s += f"    P{ip} {type} {qty_name} {values} \n"
+        return s
+
+
 class PatchHierarchy(object):
     """is a collection of patch levels"""
 
@@ -57,15 +91,14 @@ def __init__(
 
         self._sim = None
 
+        self.data_files = {}
         if data_files is not None and isinstance(data_files, dict):
             self.data_files = data_files
         elif data_files is not None:
             if hasattr(self, "data_files"):
                 self.data_files.update({data_files.filename: data_files})
             else:
                 self.data_files = {data_files.filename: data_files}
-        else:
-            self.data_files = {}
 
         self.update()
 
@@ -614,6 +647,51 @@ def dist_plot(self, **kwargs):
 
         return final, dp(final, **kwargs)
 
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        return hierarchy_array_ufunc(self, ufunc, method, *inputs, **kwargs)
+
+    def __array_function__(self, func, types, args, kwargs):
+        return hierarchy_array_function(self, func, types, args, kwargs)
+
+
+def hierarchy_array_ufunc(hier, ufunc, method, *inputs, **kwargs):
+    if method != "__call__":
+        return NotImplemented
+
+    def extract(time, ilvl):
+        return [x.level(ilvl, time) if type(x) is type(hier) else x for x in inputs]
+
+    from copy import deepcopy
+
+    ret = deepcopy(hier)
+    for time in hier.time_hier:
+        for ilvl in hier.levels(time):
+            ret.time_hier[time][ilvl] = getattr(ufunc, method)(
+                *extract(time, ilvl), **kwargs
+            )
+    return ret
+
+
+def hierarchy_array_function(hier, func, types, args, kwargs):
+    def extract(time, ilvl):
+        return [x.level(ilvl, time) if type(x) is type(hier) else x for x in args]
+
+    time_hier = {}
+    for time in hier.time_hier:
+        time_hier[time] = {}
+        for ilvl in hier.levels(time):
+            time_hier[time][ilvl] = func(*extract(time, ilvl), **kwargs)
+
+    any_level = next(iter(next(iter(time_hier.values())).values()))
+    if type(any_level[0]) is dict:
+        return time_hier  # return a dictionary of times[]
+
+    from copy import deepcopy
+
+    ret = deepcopy(hier)
+    ret.time_hier = time_hier
+    return ret
+
 
 def finest_part_data(hierarchy, time=None):
     """

pyphare/pyphare/pharesee/hierarchy/hierarchy_compute.py

@@ -0,0 +1,57 @@
+#
+#
+#
+
+import numpy as np
+
+
+def _compute_gaussian_filter_on_scalarfield(patch_datas, **kwargs):
+    from scipy.ndimage import gaussian_filter
+
+    ndim = patch_datas["value"].box.ndim
+    nb_ghosts = patch_datas["value"].ghosts_nbr[0]
+    sigma = kwargs["sigma"]
+    ds = np.asarray(patch_datas["value"][:])
+
+    ds_ = np.full(list(ds.shape), np.nan)
+
+    gf_ = gaussian_filter(ds, sigma=sigma)
+    select = tuple([slice(nb_ghosts, -nb_ghosts) for _ in range(ndim)])
+    ds_[select] = np.asarray(gf_[select])
+
+    return (
+        {"name": "value", "data": ds_, "centering": patch_datas["value"].centerings},
+    )
+
+
+def _compute_gaussian_filter_on_vectorfield(patch_datas, **kwargs):
+    from scipy.ndimage import gaussian_filter
+
+    ref_name = next(iter(patch_datas.keys()))
+
+    ndim = patch_datas[ref_name].box.ndim
+    nb_ghosts = patch_datas[ref_name].ghosts_nbr[0]
+    sigma = kwargs["sigma"]
+    ds_x = np.asarray(patch_datas["x"][:])
+    ds_y = np.asarray(patch_datas["y"][:])
+    ds_z = np.asarray(patch_datas["z"][:])
+
+    dsx_ = np.full(list(ds_x.shape), np.nan)
+    dsy_ = np.full(list(ds_y.shape), np.nan)
+    dsz_ = np.full(list(ds_z.shape), np.nan)
+
+    gfx_ = gaussian_filter(ds_x, sigma=sigma)
+    gfy_ = gaussian_filter(ds_y, sigma=sigma)
+    gfz_ = gaussian_filter(ds_z, sigma=sigma)
+
+    select = tuple([slice(nb_ghosts, -nb_ghosts) for _ in range(ndim)])
+
+    dsx_[select] = np.asarray(gfx_[select])
+    dsy_[select] = np.asarray(gfy_[select])
+    dsz_[select] = np.asarray(gfz_[select])
+
+    return (
+        {"name": "x", "data": dsx_, "centering": patch_datas["x"].centerings},
+        {"name": "y", "data": dsy_, "centering": patch_datas["y"].centerings},
+        {"name": "z", "data": dsz_, "centering": patch_datas["z"].centerings},
+    )

pyphare/pyphare/pharesee/hierarchy/hierarchy_utils.py

@@ -699,3 +699,64 @@ def _skip(qty):
                 else:
                     raise RuntimeError("unexpected state")
     return cier
+
+
+def plot_velocity_peaks_over_time(run, times, temperature, out_file, sigma=0):
+    """POC - not clean or finished"""
+    import matplotlib.pyplot as plt
+    from scipy.signal import find_peaks
+    from pyphare.pharesee.run.utils import _compute_to_primal
+
+    def get_velocities():
+        for it, t in enumerate(times):
+            Vi = run.GetVi(t)[:]
+            if sigma > 0:
+                Vi = Vi.gaussian(sigma)
+            xV = Vi.level(0, t)[0]["x"].x[2:-2]  # constant
+            v = Vi.level(0, t)[0]["x"].dataset[2:-2]  # nans on borders
+
+            if np.isnan(np.sum(v)):
+                print("summary")
+                print("nan", np.count_nonzero(np.isnan(v)))
+                print("notnan", np.count_nonzero(~np.isnan(v)))
+                raise RuntimeError("Nan!")
+
+            if it == 0:
+                vt = np.zeros((len(v), len(times)))
+            vt[:, it] = v
+
+        return xV, vt
+
+    def get_peaks(x, Vs):
+        nt = Vs.shape[1]
+        positions = np.zeros((nt, 2))
+        amplitudes = np.zeros((nt, 2))
+        for it in range(nt):
+            ps = find_peaks(Vs[:, it], height=0.010)
+            if len(ps[0]) == 1:
+                positions[it, 0] = x[ps[0][0]]
+                positions[it, 1] = x[ps[0][0]]
+                amplitudes[it, 0] = Vs[ps[0][0], it]
+                amplitudes[it, 1] = Vs[ps[0][0], it]
+            else:
+                positions[it, 0] = x[ps[0][0]]
+                positions[it, 1] = x[ps[0][1]]
+                amplitudes[it, 0] = Vs[ps[0][0], it]
+                amplitudes[it, 1] = Vs[ps[0][1], it]
+        return positions, amplitudes
+
+    x, Vs = get_velocities()
+    positions, amplitudes = get_peaks(x, Vs)
+
+    fig, ax = plt.subplots()
+
+    for it, t in enumerate(times):
+        ax.plot(x, Vs[:, it], label=r"t={:6.4f}".format(t))
+        ax.set_ylim((-0.03, 0.1))
+        ax.axhline(0, ls="--", color="k")
+        for p in positions[it, :]:
+            ax.axvline(p, color="gray", ls="-.")
+    ax.set_title("$T_e$ = {:6.4f}".format(temperature))
+    ax.legend()
+
+    fig.savefig(out_file)

pyphare/pyphare/pharesee/hierarchy/patch.py

@@ -38,7 +38,9 @@ def __repr__(self):
         return self.__str__()
 
     def __getitem__(self, key):
-        return self.patch_datas[key]
+        if key in self.patch_datas:
+            return self.patch_datas[key]
+        raise IndexError(f"No patchdata for key: {key} in {self.patch_datas}")
 
     def copy(self):
         """does not copy patchdatas.datasets (see class PatchData)"""
@@ -67,3 +69,39 @@ def __call__(self, qty, **kwargs):
             return pd.dataset[idx + nbrGhosts, nbrGhosts:-nbrGhosts]
         elif idim == 1:
             return pd.dataset[nbrGhosts:-nbrGhosts, idx + nbrGhosts]
+
+    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
+        return patch_array_ufunc(self, ufunc, method, *inputs, **kwargs)
+
+    def __array_function__(self, func, types, args, kwargs):
+        return patch_array_function(self, func, types, args, kwargs)
+
+
+def patch_array_ufunc(patch, ufunc, method, *inputs, **kwargs):
+    if method != "__call__":
+        return NotImplemented
+
+    def extract(key):
+        return [x.patch_datas[key] if type(x) is type(patch) else x for x in inputs]
+
+    patch_datas = {
+        key: getattr(ufunc, method)(*extract(key), **kwargs)
+        for key in patch.patch_datas
+    }
+    return type(patch)(
+        patch_datas, patch_id=patch.id, layout=patch.layout, attrs=patch.attrs
+    )
+
+
+def patch_array_function(patch, func, types, args, kwargs):
+    def extract(key):
+        return [x.patch_datas[key] if type(x) is type(patch) else x for x in args]
+
+    final = {key: func(*extract(key), **kwargs) for key in patch.patch_datas}
+    any_patch_data = next(iter(patch.patch_datas.values()))
+    any_data = next(iter(final.values()))
+
+    if type(any_patch_data) is not type(any_data):
+        return final
+
+    return type(patch)(final, patch_id=patch.id, layout=patch.layout, attrs=patch.attrs)