From 6e182ebb358f81ed8208752e17914170ff0dfb62 Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Tue, 16 Mar 2021 17:06:29 +0100 Subject: [PATCH] [postprocessor] [plotter] Make generic rules averageables + cleanup --- plotter.py | 147 ++++++++++++------------------ postprocessor.py | 226 +++++++++++++++++++++-------------------------- 2 files changed, 160 insertions(+), 213 deletions(-) diff --git a/plotter.py b/plotter.py index a639e41..83af97e 100644 --- a/plotter.py +++ b/plotter.py @@ -1269,39 +1269,6 @@ class Plotter(Aggregator, BaseProcessor): "Map of the grav Shakura&Sunaev alpha parameter for disks", dependencies=["alpha_grav"], ), - "vphi": PlotRule( - self, - partial( - self._plot_map, - "vphi", - label=r"$v_\phi$", - # unit=U.km_s - ), - "Azimuthal speed", - dependencies=["vphi"], - ), - "vr": PlotRule( - self, - partial( - self._plot_map, - "vr", - label=r"$v_r$", - # unit=U.km_s - ), - "Radial speed", - dependencies=["vr"], - ), - "rho": PlotRule( - self, - partial( - self._plot_map, - "rho", - label=r"$\rho$", - # unit=U.Msun_pc3 - ), - "Density slice at s = 0, with s = x, y or z.", - dependencies=["rho"], - ), "coldens_l": PlotRule( self, partial( @@ -1314,41 +1281,41 @@ class Plotter(Aggregator, BaseProcessor): "Column density with level overlay", dependencies=["coldens", "levels"], ), - "rho_v": PlotRule( + "slice_rho_v": PlotRule( self, partial( self._plot_map, - "rho", + "slice_rho", label=r"$\rho$", unit=U.Msun_pc3, overlays=[self._overlay_speed], ), "Density slice with speed overlay", - dependencies=["rho", "speed_h", "speed_v"], + dependencies=["slice_rho", "speed_h", "speed_v"], ), - "rho_B": PlotRule( + "slice_rho_B": PlotRule( self, partial( self._plot_map, - "rho", + "slice_rho", label=r"$\rho$", unit=U.Msun_pc3, overlays=[self._overlay_B], ), "Density slice with magnetic field overlay", - dependencies=["rho", "B_h", "B_v"], + dependencies=["slice_rho", "B_h", "B_v"], ), - "rho_B_vel": PlotRule( + "slice_rho_B_vel": PlotRule( self, partial( self._plot_map, - "rho", + "slice_rho", label=r"$\rho$", unit=U.Msun_pc3, overlays=[self._overlay_B, self._overlay_speed], ), "Density slice with magnetic field and velocity overlay", - dependencies=["rho", "B_h", "B_v", "speed_h", "speed_v"], + dependencies=["slice_rho", "B_h", "B_v", "speed_h", "speed_v"], ), "jeans_ratio": PlotRule( self, @@ -1629,17 +1596,58 @@ class Plotter(Aggregator, BaseProcessor): averageables = [ "coldens", - "rho", - "T", "Q", - "vr", - "vphi", - "rho_avg", - "P_avg", - "T_avg", + "T", + "T_mwavg", "alpha_disk", "alpha_grav", ] + + # Generic rules directly from Ramses fields + for field in self.pp_params.pymses.variables: + + def generic_rule(name): + + self.rules["slice_" + name] = PlotRule( + self, + partial(self._plot_map, "slice_" + name), + "{} slice".format(name), + dependencies=["slice_" + name], + ) + + self.rules[name + "_mwavg"] = PlotRule( + self, + partial(self._plot_map, name + "_mwavg"), + "Ax mass-weighted averaged {}".format(name), + dependencies=[name + "_mwavg"], + ) + + self.rules[name + "_avg"] = PlotRule( + self, + partial(self._plot_map, name + "_avg"), + "Ax averaged {}".format(name), + dependencies=[name + "_avg"], + ) + + averageables.append("slice_" + name) + averageables.append(name + "_mwavg") + averageables.append(name + "_avg") + + # special for vectors + if field in ["g", "vel"]: + # Components + for i, dir in enumerate(["x", "y", "z"]): + generic_rule(field + dir) + + # Radial + generic_rule(field + "r") + # Orthoradial + generic_rule(field + "phi") + # Norm + generic_rule(field + "_norm") + else: + generic_rule(field) + for name in averageables: self.rules["rad_" + name] = PlotRule( self, @@ -1698,47 +1706,6 @@ class Plotter(Aggregator, BaseProcessor): for name in ["time", "dt", "a", "mempc1", "mempc2"]: self._gen_from_log("fine_step_from_log", name_y=name, name_x="fine_step") - # Generic rules directly from Ramses fields - for field in self.pp_params.pymses.variables: - - def generic_rule(name): - - self.rules["slice_" + name] = PlotRule( - self, - partial(self._plot_map, "slice_" + name), - "{} slice".format(name), - dependencies=["slice_" + name], - ) - - self.rules[name + "_mwavg"] = PlotRule( - self, - partial(self._plot_map, name + "_mwavg"), - "Ax mass-weighted averaged {}".format(name), - dependencies=[name + "_mwavg"], - ) - - self.rules[name + "_avg"] = PlotRule( - self, - partial(self._plot_map, name + "_avg"), - "Ax averaged {}".format(name), - dependencies=[name + "_avg"], - ) - - # special for vectors - if field in ["g", "vel"]: - # Components - for i, dir in enumerate(["x", "y", "z"]): - generic_rule(field + dir) - - # Radial - generic_rule(field + "r") - # Orthoradial - generic_rule(field + "phi") - # Norm - generic_rule(field + "_norm") - else: - generic_rule(field) - # Dict of overlays self.overlays = { "B": self._overlay_B, diff --git a/postprocessor.py b/postprocessor.py index 1848111..ad3cddc 100644 --- a/postprocessor.py +++ b/postprocessor.py @@ -770,10 +770,6 @@ class PostProcessor(HDF5Container): datamap = self._rt.process(self._cam[ax_los], surf_qty=True) return datamap.map.T - def _rho(self, ax_los, z=0.0): - datamap_rho = slicing.SliceMap(self._amr, self._cam[ax_los], self._rho_op, z=z) - return (datamap_rho.map).T - def _vector_h(self, name, unit, ax_los, z=0.0): h_op = ScalarOperator( lambda dset: dset[name][:, self._ax_nb[self._axes_h[ax_los]]], @@ -816,7 +812,7 @@ class PostProcessor(HDF5Container): def _temperature(self, ax_los, z=0.0): P_op = ScalarOperator(lambda dset: dset["P"], self._ro.info["unit_pressure"]) dmap_P = (slicing.SliceMap(self._amr, self._cam[ax_los], P_op, z=z)).map.T - dmap_rho = self.save.get_node("/maps/rho_{}".format(ax_los)).read() + dmap_rho = self.save.get_node("/maps/slice_rho_{}".format(ax_los)).read() return dmap_P / dmap_rho def _levels(self, ax_los): @@ -830,7 +826,7 @@ class PostProcessor(HDF5Container): def _jeans(self, ax_los): dmap_T = self.save.get_node("/maps/T_" + ax_los).read() - dmap_rho = self.save.get_node("/maps/rho_" + ax_los).read() + dmap_rho = self.save.get_node("/maps/slice_rho_" + ax_los).read() dmap_jeans = np.sqrt(np.pi * dmap_T / dmap_rho) return dmap_jeans @@ -1108,10 +1104,10 @@ class PostProcessor(HDF5Container): radial_bins = self.save.get_node("/radial/radial_bins_" + ax_los).read() mean_bin_vr = self.save.get_node( - "/radial/rad_avg_" + "vr" + "_" + ax_los + "/radial/rad_avg_" + "slice_velr" + "_" + ax_los ).read() mean_bin_vphi = self.save.get_node( - "/radial/rad_avg_" + "vphi" + "_" + ax_los + "/radial/rad_avg_" + "slice_velphi" + "_" + ax_los ).read() mean_bin_vr = np.concatenate((mean_bin_vr, [mean_bin_vr[-1]])) mean_bin_vphi = np.concatenate((mean_bin_vphi, [mean_bin_vr[-1]])) @@ -1324,7 +1320,7 @@ class PostProcessor(HDF5Container): # Equilibrium Gvry, Gvrx = np.gradient(vr, dy, dx, edge_order=2) gradvr = (xx * Gvrx + yy * Gvry) / rr - dvr = gradvr + vr * gradvr # Complete derivative + dvr = vr * gradvr # Complete derivative # Thermal support GPy, GPx = np.gradient(Pz, dy, dx, edge_order=2) @@ -1363,30 +1359,6 @@ class PostProcessor(HDF5Container): "/maps", unit=self.info["unit_density"] * self.info["unit_length"], ), - "vr": Rule( - self, - partial( - self._ax_avg, - self.oct_getter_vr, - mass_weighted=True, - unit=self.info["unit_velocity"], - ), - "Vertically mass-weighted averaged radial velocity", - "/maps", - unit=self.info["unit_velocity"], - ), - "vphi": Rule( - self, - partial( - self._ax_avg, - self.oct_getter_vphi, - mass_weighted=True, - unit=self.info["unit_velocity"], - ), - "Vertically mass-weighted averaged azimuthal velocity", - "/maps", - unit=self.info["unit_velocity"], - ), "T_mwavg": Rule( self, partial( @@ -1414,10 +1386,10 @@ class PostProcessor(HDF5Container): "/maps", unit=U.none, dependencies=[ - "avg_map_rho_avg", + "avg_map_slice_rho_avg", "avg_map_T_mwavg", - "avg_map_vr", - "avg_map_vphi", + "avg_map_slice_velr", + "avg_map_slice_velphi", ], ), "alpha_grav": Rule( @@ -1429,13 +1401,6 @@ class PostProcessor(HDF5Container): unit=U.none, dependencies=["avg_map_coldens", "avg_map_T_mwavg"], ), - "rho": Rule( - self, - self._rho, - "Density slice", - "/maps", - unit=self.info["unit_density"], - ), "speed_h": Rule( self, self._speed_h, @@ -1452,14 +1417,14 @@ class PostProcessor(HDF5Container): ), "B_h": Rule( self, - self._speed_h, + self._B_h, "Horizontal slice of the magnetic field wrt the line of sight", "/maps", unit=self.info["unit_mag"], ), "B_v": Rule( self, - self._speed_v, + self._B_v, "Vertical slice of the magnetic field wrt the line of sight", "/maps", unit=self.info["unit_mag"], @@ -1469,7 +1434,7 @@ class PostProcessor(HDF5Container): self._temperature, "Temperature slice", "/maps", - dependencies=["rho"], + dependencies=["slice_rho"], unit=self.info["unit_pressure"] / self.info["unit_density"], ), "levels": Rule( @@ -1480,7 +1445,7 @@ class PostProcessor(HDF5Container): self._jeans, "Jeans length slice", "/maps", - dependencies=["rho", "T"], + dependencies=["slice_rho", "T"], ), "jeans_ratio": Rule( self, @@ -1550,6 +1515,7 @@ class PostProcessor(HDF5Container): self._radial_bins, "Radial bins", "/radial", + unit=self.info["unit_length"] / self.lbox, ), "radial_centers": Rule( self, @@ -1557,8 +1523,15 @@ class PostProcessor(HDF5Container): "Centers of radial bins", "/radial", dependencies=["radial_bins"], + unit=self.info["unit_length"] / self.lbox, + ), + "rr": Rule( + self, + self._rr, + "Coordinate map", + "/maps", + unit=self.info["unit_length"] / self.lbox, ), - "rr": Rule(self, self._rr, "Coordinate map", "/maps"), "bins_on_map": Rule( self, self._bins_on_map, @@ -1571,7 +1544,7 @@ class PostProcessor(HDF5Container): self._B_int, "Magnetic intensity slice", "/maps", - dependencies=["rho"], + dependencies=["slice_rho"], unit=self.info["unit_mag"], ), # PDF @@ -1686,29 +1659,96 @@ class PostProcessor(HDF5Container): ), } - # Average and other averageables = [ "coldens", - "rho", - "T", "Q", - "vphi", - "vr", - "rho_avg", - "P_avg", - "T_avg", - "P_mwavg", + "T", "T_mwavg", "alpha_disk", "alpha_grav", ] + + # Generic rules directly from Ramses fields + for field in self.pp_params.pymses.variables: + + def generic_rule(name, getter, unit, oct_getter=None): + if oct_getter is None: + oct_getter = getter + + self.rules["slice_" + name] = Rule( + self, + partial(self._slice, getter, z=0.0, unit=unit), + "{} slice".format(name), + "/maps", + unit=unit, + ) + + self.rules[name + "_mwavg"] = Rule( + self, + partial(self._ax_avg, oct_getter, mass_weighted=True, unit=unit), + "Ax mass-weighted averaged {}".format(name), + "/maps", + unit=unit, + ) + + self.rules[name + "_avg"] = Rule( + self, + partial(self._ax_avg, oct_getter, mass_weighted=False, unit=unit), + "Ax averaged {}".format(name), + "/maps", + unit=unit, + ) + + averageables.append("slice_" + name) + averageables.append(name + "_mwavg") + averageables.append(name + "_avg") + + # special for vectors + if field in ["g", "vel"]: + # Components + for i, dir in enumerate(["x", "y", "z"]): + generic_rule( + field + dir, + partial(vect_getter, field, i), + self.unit_key[field], + ) + + # Radial + generic_rule( + field + "r", + partial(self.getter_vect_r, name_vect=field), + self.unit_key[field], + oct_getter=partial(self.oct_getter_vect_r, name_vect=field), + ) + + # Othoradial + generic_rule( + field + "phi", + partial(self.getter_vect_phi, name_vect=field), + self.unit_key[field], + oct_getter=partial(self.oct_getter_vect_phi, name_vect=field), + ) + + # Norm + generic_rule( + field + "_norm", partial(norm_getter, field), self.unit_key[field] + ) + + else: + generic_rule(field, partial(simple_getter, field), self.unit_key[field]) + + # radial average and other for name in averageables: + + unit = self.rules[name].unit + self.rules["rad_avg_" + name] = Rule( self, partial(self._rad_avg, name), "Azimuthal average of {}".format(name), "/radial", dependencies=["radial_centers", "bins_on_map", name], + unit=unit, ) self.rules["rad_mwavg_" + name] = Rule( self, @@ -1716,6 +1756,7 @@ class PostProcessor(HDF5Container): "Mass weighted azimuthal average of {}".format(name), "/radial", dependencies=["coldens", "radial_centers", "bins_on_map", name], + unit=unit, ) self.rules["avg_map_" + name] = Rule( self, @@ -1723,6 +1764,7 @@ class PostProcessor(HDF5Container): "Interpolated map of azimuthal average of {}".format(name), "/maps", dependencies=["radial_centers", "bins_on_map", "rad_avg_" + name], + unit=unit, ) self.rules["mwavg_map_" + name] = Rule( self, @@ -1730,6 +1772,7 @@ class PostProcessor(HDF5Container): "Interpolated map of azimuthal average of {}".format(name), "/maps", dependencies=["radial_centers", "bins_on_map", "rad_mwavg_" + name], + unit=unit, ) self.rules["fluct_" + name] = Rule( self, @@ -1769,6 +1812,7 @@ class PostProcessor(HDF5Container): dependencies=["pdf_" + name], ) for name_bin in averageables: + unit_bin = self.rules[name_bin].unit if name_bin is not name: self.rules["mbb_" + name + "_" + name_bin] = Rule( self, @@ -1776,75 +1820,11 @@ class PostProcessor(HDF5Container): "Mean of {} by bins of {}".format(name, name_bin), "/hist", dependencies=[name, name_bin], + unit=[unit, unit_bin], ) self._gen_rule_transform("fluct_coldens", np.nanmax, "max", group="/globals") - # Generic rules directly from Ramses fields - for field in self.pp_params.pymses.variables: - - def generic_rule(name, getter, unit, oct_getter=None): - if oct_getter is None: - oct_getter = getter - - self.rules["slice_" + name] = Rule( - self, - partial(self._slice, getter, z=0.0, unit=unit), - "{} slice".format(name), - "/maps", - unit=unit, - ) - - self.rules[name + "_mwavg"] = Rule( - self, - partial(self._ax_avg, oct_getter, mass_weighted=True, unit=unit), - "Ax mass-weighted averaged {}".format(name), - "/maps", - unit=unit, - ) - - self.rules[name + "_avg"] = Rule( - self, - partial(self._ax_avg, oct_getter, mass_weighted=False, unit=unit), - "Ax averaged {}".format(name), - "/maps", - unit=unit, - ) - - # special for vectors - if field in ["g", "vel"]: - # Components - for i, dir in enumerate(["x", "y", "z"]): - generic_rule( - field + dir, - partial(vect_getter, field, i), - self.unit_key[field], - ) - - # Radial - generic_rule( - field + "r", - partial(self.getter_vect_r, name_vect=field), - self.unit_key[field], - oct_getter=self.oct_getter_vect_r, - ) - - # Othoradial - generic_rule( - field + "phi", - partial(self.getter_vect_phi, name_vect=field), - self.unit_key[field], - oct_getter=self.oct_getter_vect_phi, - ) - - # Norm - generic_rule( - field + "_norm", partial(norm_getter, field), self.unit_key[field] - ) - - else: - generic_rule(field, partial(simple_getter, field), self.unit_key[field]) - super(PostProcessor, self).def_rules()