[postprocessor] [plotter] Make generic rules averageables + cleanup

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",
+            "T",
-            "vphi",
+            "T_mwavg",
-            "rho_avg",
-            "P_avg",
-            "T_avg",
             "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,

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_slice_velr",
-                    "avg_map_vphi",
+                    "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",
+            "T",
-            "vr",
-            "rho_avg",
-            "P_avg",
-            "T_avg",
-            "P_mwavg",
             "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()