From 5216e3953a7de92dd38153597318e7fd877cdd24 Mon Sep 17 00:00:00 2001 From: Nicolas Aunai Date: Mon, 4 Oct 2021 13:46:28 +0200 Subject: [PATCH] improve plotting (#586) current 2D, plot2d patches, yee finest coord --- pyphare/pyphare/core/gridlayout.py | 44 +++++++++---- pyphare/pyphare/pharesee/hierarchy.py | 92 +++++++++++++++++++++++++-- pyphare/pyphare/pharesee/plotting.py | 28 +++++++- pyphare/pyphare/pharesee/run.py | 86 +++++++++++++++++++++---- 4 files changed, 215 insertions(+), 35 deletions(-) diff --git a/pyphare/pyphare/core/gridlayout.py b/pyphare/pyphare/core/gridlayout.py index 857fdbab5..dc17663e0 100644 --- a/pyphare/pyphare/core/gridlayout.py +++ b/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 diff --git a/pyphare/pyphare/pharesee/hierarchy.py b/pyphare/pyphare/pharesee/hierarchy.py index cfaf1dcb4..9e93f1325 100644 --- a/pyphare/pyphare/pharesee/hierarchy.py +++ b/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) diff --git a/pyphare/pyphare/pharesee/plotting.py b/pyphare/pyphare/pharesee/plotting.py index a7ee6cae9..d1957f629 100644 --- a/pyphare/pyphare/pharesee/plotting.py +++ b/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, diff --git a/pyphare/pyphare/pharesee/run.py b/pyphare/pyphare/pharesee/run.py index 5185a778a..d8b3114b9 100644 --- a/pyphare/pyphare/pharesee/run.py +++ b/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