improve plotting (#586)

current 2D, plot2d patches, yee finest coord

pyphare/pyphare/core/gridlayout.py

@@ -58,6 +58,34 @@ def __init__(self):
         self.centerY = yee_centering["y"]
         self.centerZ = yee_centering["z"]
 
+
+def yeeCoordsFor(origin, nbrGhosts, dl, nbrCells, qty, direction, withGhosts=False):
+
+    assert direction in direction_to_dim, f"direction ({direction} not supported)"
+    assert qty in yee_centering[direction] or qty in yee_centering_lower[direction], f"qty ({qty} not supported)"
+    if qty in yee_centering_lower[direction] and qty not in yee_centering[direction]:
+        qty = qty[0].upper() + qty[1:]
+
+    centering = yee_centering[direction][qty]
+
+    offset = 0
+    dim = direction_to_dim[direction]
+    if withGhosts:
+        size = nbrCells[dim] + (nbrGhosts * 2)
+    else:
+        size = nbrCells[dim]
+
+    if centering == 'dual':
+        offset = 0.5*dl[dim]
+    else:
+        size += 1
+
+    if withGhosts:
+        return origin[dim] - nbrGhosts * dl[dim] + np.arange(size) * dl[dim] + offset
+    else:
+        return origin[dim] + np.arange(size) * dl[dim] + offset
+
+
 class GridLayout(object):
 
     def __init__(self, box=Box(0,0), origin=0, dl=0.1, interp_order=1):
@@ -256,24 +284,12 @@ def yeeCoordsFor(self, qty, direction):
         :param direction: can only be a single one
         """
 
-        assert direction in direction_to_dim, f"direction ({direction} not supported)"
-        assert qty in yee_centering[direction] or qty in yee_centering_lower[direction], f"qty ({qty} not supported)"
         if qty in yee_centering_lower[direction] and qty not in yee_centering[direction]:
             qty = qty[0].upper() + qty[1:]
 
-        centering = yee_centering[direction][qty]
-        nbrGhosts = self.nbrGhosts(self.interp_order, centering)
-
-        offset = 0
-        dim = direction_to_dim[direction]
-        size = self.box.shape[dim] + (nbrGhosts * 2)
-
-        if centering == 'dual':
-            offset = 0.5*self.dl[dim]
-        else:
-            size += 1
+        return yeeCoordsFor(self.origin, self.nbrGhosts(self.interp_order, yee_centering[direction][qty]),
+                            self.dl, self.box.shape, qty, direction, withGhosts=True)
 
-        return self.origin[dim] - nbrGhosts * self.dl[dim] + np.arange(size) * self.dl[dim] + offset
 
 
 

pyphare/pyphare/pharesee/hierarchy.py

@@ -266,7 +266,7 @@ def overlap_mask_1d(x, dl, level, qty):
             overlaped_idx = np.where( (x > xmin) & (x < xmax) )[0]
 
             is_overlaped[overlaped_idx] = True
- 
+
         else:
             raise ValueError("level needs to have finer grid resolution than that of x")
 
@@ -541,10 +541,14 @@ def level(self, level_number, time=None):
         return self.levels(time)[level_number]
 
 
-    def levelNbr(self, time):
+    def levelNbr(self, time=None):
+        if time is None:
+            time = self._default_time()
         return len(self.levels(time).items())
 
-    def levelNbrs(self, time):
+    def levelNbrs(self, time=None):
+        if time is None:
+            time = self._default_time()
         return list(self.levels(time).keys())
 
     def is_homogeneous(self):
@@ -659,7 +663,7 @@ def plot_2d_patches(self, ilvl, collections, **kwargs):
 
         return fig
 
-    def plot(self, **kwargs):
+    def plot1d(self, **kwargs):
         """
         plot
         """
@@ -705,6 +709,81 @@ def plot(self, **kwargs):
         if "filename" in kwargs:
             fig.savefig(kwargs["filename"])
 
+
+    def plot2d(self, **kwargs):
+        from mpl_toolkits.axes_grid1 import make_axes_locatable
+        from matplotlib.patches import Rectangle
+        time = kwargs.get("time", self._default_time())
+        usr_lvls = kwargs.get("levels",self.levelNbrs(time))
+        qty = kwargs.get("qty", None)
+
+        if "ax" not in kwargs:
+            fig, ax = plt.subplots()
+        else:
+            ax = kwargs["ax"]
+            fig = ax.figure
+
+        for lvl_nbr, lvl  in self.levels(time).items():
+            if lvl_nbr not in usr_lvls:
+                continue
+            for patch in self.level(lvl_nbr, time).patches:
+                pdat = patch.patch_datas[qty]
+                data = pdat.dataset[:]
+                nbrGhosts = pdat.ghosts_nbr
+                x = pdat.x
+                y = pdat.y
+                nx,ny =x.size, y.size
+                data = data.reshape((nx,ny))
+                data = data[nbrGhosts[0]:-nbrGhosts[0], nbrGhosts[1]:-nbrGhosts[1]]
+                x = x[nbrGhosts[0]:-nbrGhosts[0]]
+                y = y[nbrGhosts[1]:-nbrGhosts[1]]
+                dx,dy = pdat.layout.dl
+                x -= dx*0.5
+                y -= dy*0.5
+                x = np.append(x, x[-1]+dx)
+                y = np.append(y, y[-1]+dy)
+                im = ax.pcolormesh(x, y, data.T,
+                                   cmap=kwargs.get("cmap","Spectral_r"),
+                                   vmin=kwargs.get("vmin", None),
+                                   vmax=kwargs.get("vmax", None))
+
+                if kwargs.get("plot_patches", False) is True:
+                    r = Rectangle((patch.box.lower[0]*dx,
+                                   patch.box.lower[1]*dy),
+                                  patch.box.shape[0]*dx,
+                                  patch.box.shape[1]*dy,
+                                  fc="none", ec="k",
+                                  alpha=0.4, lw=0.8)
+                    ax.add_patch(r)
+
+        ax.set_aspect("equal")
+        ax.set_title(kwargs.get("title", ""))
+        ax.set_xlabel(kwargs.get("xlabel", "x"))
+        ax.set_ylabel(kwargs.get("ylabel", "y"))
+        if "xlim" in kwargs:
+            ax.set_xlim(kwargs["xlim"])
+        if "ylim" in kwargs:
+            ax.set_ylim(kwargs["ylim"])
+
+        divider = make_axes_locatable(ax)
+        cax = divider.append_axes("right", size="5%", pad=0.08)
+        plt.colorbar( im, ax=ax, cax=cax )
+
+        if kwargs.get("legend", None) is not None:
+            ax.legend()
+
+        if "filename" in kwargs:
+            fig.savefig(kwargs["filename"])
+
+        return fig,ax
+
+    def plot(self, **kwargs):
+        if self.ndim == 1:
+            return self.plot1d(**kwargs)
+        elif self.ndim==2:
+            return self.plot2d(**kwargs)
+
+
     def dist_plot(self, **kwargs):
         """
         plot phase space of a particle hierarchy
@@ -927,6 +1006,7 @@ def compute_hier_from(h, compute):
 
      caveat: routine only works in 1D so far.
     """
+    assert(len(h.time_hier) == 1) # only single time hierarchies now
     patch_levels = {}
     for ilvl, lvl in h.patch_levels.items():
         patches = {}
@@ -945,7 +1025,9 @@ def compute_hier_from(h, compute):
 
         patch_levels[ilvl] = PatchLevel(ilvl, patches[ilvl])
 
-    return PatchHierarchy(patch_levels, h.domain_box, refinement_ratio)
+    t = list(h.time_hier.keys())[0]
+    return PatchHierarchy(patch_levels, h.domain_box, refinement_ratio,
+                          time=t)
 
 
 

pyphare/pyphare/pharesee/plotting.py

@@ -218,6 +218,7 @@ def finest_field_plot(run_path, qty, **kwargs):
     from pyphare.pharesee.hierarchy import get_times_from_h5
     from pyphare.pharesee.run import Run
     from mpl_toolkits.axes_grid1 import make_axes_locatable
+    import pyphare.core.gridlayout as gridlayout
 
     r = Run(run_path)
 
@@ -254,6 +255,13 @@ def finest_field_plot(run_path, qty, **kwargs):
             time = times[0]
         interpolator, finest_coords = r.GetNi(time, merged=True,\
                                               interp=interp)[qty]
+    elif qty in ("Jx", "Jy", "Jz"):
+        file = os.path.join(run_path, "EM_B.h5")
+        if time is None:
+            times = get_times_from_h5(file)
+            time = times[0]
+        interpolator, finest_coords = r.GetJ(time, merged=True,\
+                                              interp=interp)[qty]
     else:
         # ___ TODO : should also include the files for a given population
         raise ValueError("qty should be in ['Bx', 'By', 'Bz', 'Ex', 'Ey', 'Ez', 'Fx', 'Fy', 'Fz', 'Vx', 'Vy', 'Vz', 'rho']")
@@ -270,14 +278,28 @@ def finest_field_plot(run_path, qty, **kwargs):
         ax.plot(finest_coords[0], interpolator(finest_coords[0]),\
                 drawstyle=drawstyle)
     elif dim == 2:
-        X, Y = np.meshgrid(finest_coords[0], finest_coords[1])
+
+
+        x = finest_coords[0]
+        y = finest_coords[1]
+        dx = x[1] - x[0]
+        dy = y[1] - y[0]
+
+        # pcolormesh considers DATA_ij to be the center of the pixel
+        # and X,Y are the corners so XY need to be made 1 value larger
+        # and shifted around DATA_ij
+        x -= dx/2
+        x= np.append(x, x[-1]+dx)
+        y -= dy/2
+        y= np.append(y, y[-1]+dy)
+
+        X, Y = np.meshgrid(x, y)
         DATA = interpolator(X, Y)
 
         vmin = kwargs.get("vmin", np.nanmin(DATA))
         vmax = kwargs.get("vmax", np.nanmax(DATA))
         cmap = kwargs.get("cmap", 'Spectral_r')
-
-        im = ax.pcolormesh(finest_coords[0], finest_coords[1],
+        im = ax.pcolormesh(x,y,
                    DATA,
                    cmap = cmap,
                    vmin = vmin,

pyphare/pyphare/pharesee/run.py

@@ -24,18 +24,72 @@ def _current1d(by, bz, xby, xbz):
     jz[-1]=jz[-2]
     return jy, jz
 
+def _current2d(bx, by, bz, dx, dy):
+    # jx = dyBz
+    # jy = -dxBz
+    # jz = dxBy - dyBx
+    # the following hard-codes yee layout
+    # which is not true in general
+    # we should at some point provide proper
+    # derivation routines in the gridlayout
+    jx = np.zeros(by.shape)
+    jy = np.zeros(bx.shape)
+    jz = np.zeros((bx.shape[0], by.shape[1]))
+
+    jx[:,1:-1] = (bz[:,1:] - bz[:,:-1])/dy
+    jy[1:-1,:] = -(bz[1:,:] - bz[:-1,:])/dx
+    jz[1:-1,1:-1] = (by[1:, 1:-1] - by[:-1 , 1:-1])/dx - (bx[1:-1 , 1:] - bx[1:-1, :-1])/dy
+
+    jy[0,:]  = jy[1,:]
+    jy[:,0]  = jy[:,1]
+    jy[-1,:] = jy[-2,:]
+    jy[:,-1] = jy[:,-2]
+
+    jz[0,:]  = jz[1,:]
+    jz[:,0]  = jz[:,1]
+    jz[-1,:] = jz[-2,:]
+    jz[:,-1] = jz[:,-2]
+
+    jx[0,:]  = jx[1,:]
+    jx[:,0]  = jx[:,1]
+    jx[-1,:] = jx[-2,:]
+    jx[:,-1] = jx[:,-2]
+
+    return jx, jy, jz
+
+
 
 def _compute_current(patch):
-    By = patch.patch_datas["By"].dataset[:]
-    xby  = patch.patch_datas["By"].x
-    Bz = patch.patch_datas["Bz"].dataset[:]
-    xbz  = patch.patch_datas["Bz"].x
-    Jy, Jz =  _current1d(By, Bz, xby, xbz)
-    return ({"name":"Jy", "data":Jy,"centering":"primal"},
-            {"name":"Jz", "data":Jz,"centering":"primal"})
+    if patch.box.ndim ==1 : 
+        By = patch.patch_datas["By"].dataset[:]
+        xby  = patch.patch_datas["By"].x
+        Bz = patch.patch_datas["Bz"].dataset[:]
+        xbz  = patch.patch_datas["Bz"].x
+        Jy, Jz =  _current1d(By, Bz, xby, xbz)
+        return ({"name":"Jy", "data":Jy,"centering":"primal"},
+                {"name":"Jz", "data":Jz,"centering":"primal"})
+
+    elif patch.box.ndim ==2 :
+        Bx = patch.patch_datas["Bx"].dataset[:]
+        By = patch.patch_datas["By"].dataset[:]
+        Bz = patch.patch_datas["Bz"].dataset[:]
+
+        dx,dy  = patch.dl
 
+        Jx, Jy, Jz =  _current2d(Bx, By, Bz, dx, dy)
+        #return ({"name":"Jx", "data":Jx,"centering":"primal"},
+        #        {"name":"Jy", "data":Jy,"centering":"primal"},
+        #        {"name":"Jz", "data":Jz,"centering":"primal"})
 
-def make_interpolator(data, coords, interp, domain, dl):
+
+        components = ("Jx", "Jy", "Jz")
+        centering = {component:[patch.layout.centering[direction][component] for direction in ("X", "Y")] for component in components}
+
+        return ({"name":"Jx", "data":Jx, "centering":centering["Jx"]},
+                {"name":"Jy", "data":Jy, "centering":centering["Jy"]},
+                {"name":"Jz", "data":Jz,"centering":centering["Jz"]})
+
+def make_interpolator(data, coords, interp, domain, dl, qty, nbrGhosts):
     """
     :param data: the values of the data that will be used for making
     the interpolator, defined on coords
@@ -45,7 +99,7 @@ def make_interpolator(data, coords, interp, domain, dl):
     finest_coords will be the structured coordinates defined on the
     finest grid.
     """
-
+    from pyphare.core.gridlayout import yeeCoordsFor
     dim = coords.ndim
 
     if dim == 1:
@@ -57,8 +111,8 @@ def make_interpolator(data, coords, interp, domain, dl):
                                 assume_sorted=False)
 
         nx = 1+int(domain[0]/dl[0])
-        x = dl[0]*np.arange(0, nx)
 
+        x = yeeCoordsFor([0]*dim, nbrGhosts, dl, [nx],  qty, "x")
         finest_coords = (x,)
 
     elif dim == 2:
@@ -72,8 +126,11 @@ def make_interpolator(data, coords, interp, domain, dl):
         else:
             raise ValueError("interp can only be 'nearest' or 'bilinear'")
 
-        x = np.arange(0, domain[0]+dl[0], dl[0])
-        y = np.arange(0, domain[1]+dl[1], dl[1])
+        nCells = [1+int(d/dl) for d,dl in zip(domain, dl)]
+        x = yeeCoordsFor([0]*dim, nbrGhosts, dl, nCells,  qty, "x")
+        y = yeeCoordsFor([0]*dim, nbrGhosts, dl, nCells,  qty, "y")
+        #x = np.arange(0, domain[0]+dl[0], dl[0])
+        #y = np.arange(0, domain[1]+dl[1], dl[1])
         finest_coords = (x, y)
 
     else:
@@ -101,11 +158,14 @@ def _get(self, hierarchy, time, merged, interp):
             domain = self.GetDomainSize()
             dl = self.GetDl()
 
+            # assumes all qties in the hierarchy have the same ghost width
+            # so take the first patch data of the first patch of the first level....
+            nbrGhosts = list(hierarchy.level(0).patches[0].patch_datas.values())[0].ghosts_nbr
             merged_qties = {}
             for qty in hierarchy.quantities():
                 data, coords = flat_finest_field(hierarchy, qty, time=time)
                 merged_qties[qty] = make_interpolator(data, coords,\
-                                                      interp, domain, dl)
+                                                      interp, domain, dl, qty, nbrGhosts)
             return merged_qties
         else:
             return hierarchy