From d92fb888f9c8c0c52f79eb26bc7bf96546e5a765 Mon Sep 17 00:00:00 2001 From: Hennebelle Patrick Date: Mon, 22 Jun 2020 14:57:08 +0200 Subject: [PATCH] Magnetic Intensity map, streamlines and averaged field versus density --- plotter.py | 32 ++++++++++++++++++++++ postprocessor.py | 70 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 102 insertions(+) diff --git a/plotter.py b/plotter.py index aef79e2..bc8e54c 100644 --- a/plotter.py +++ b/plotter.py @@ -987,6 +987,18 @@ class Plotter(Aggregator, BaseProcessor): "Density slice with magnetic field overlay", dependencies=["rho", "B_h", "B_v"], ), + "rho_B_vel": PlotRule( + self, + partial( + self._plot_map, + "rho", + label=r"$\rho$", + unit=cst.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"], + ), "jeans_ratio": PlotRule( self, partial( @@ -1019,6 +1031,12 @@ class Plotter(Aggregator, BaseProcessor): "$\rho$-PDF", dependencies=["rho_pdf"], ), + "cos_pdf": PlotRule( + self, + partial(self._plot_hist, "cos_pdf"), + "cos-PDF", + dependencies=["cos_pdf"], + ), "avg_coldens_pdf": PlotRule( self, partial( @@ -1044,6 +1062,20 @@ class Plotter(Aggregator, BaseProcessor): "P-PDF on a 2D slice ", dependencies=["P_pdf"], ), + "B_int": PlotRule( + self, + partial( + self._plot_map, "B_int", label=r"$\mid \mathrm{B} \mid$", unit=cst.T + ), + "Magnetic intensity map", + dependencies=["B_int"], + ), + "Brho": PlotRule( + self, + partial(self._plot, "/datasets/Brho/rho", "/datasets/Brho/B"), + "Brho on a 2D slice ", + dependencies=["Brho"], + ), "rho_prof": PlotRule( self, partial(self._plot, "/profile/axis", "/profile/rho_prof"), diff --git a/postprocessor.py b/postprocessor.py index f36198d..1b556cb 100644 --- a/postprocessor.py +++ b/postprocessor.py @@ -8,6 +8,23 @@ vol_func = lambda dset: dset["dx"] ** 3 # Volume function getter_T = lambda dset: dset["P"] / dset["rho"] # Temperature +def _getter_abs_cos_vB(dset): + B_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1)) + v_norm = np.sqrt(np.sum(dset["vel"] ** 2, axis=1)) + # Compute the dot product in each cell + dot_prod = np.einsum("ij,ij->i", dset["vel"], dset["Br"]) + return np.abs(dot_prod) / (v_norm * B_norm) + + +def _getter_B_int(dset): + B_norm = np.sqrt(np.sum(dset["Br"] ** 2, axis=1)) + return B_norm + + +def _getter_rho(dset): + return dset["rho"] + + class PostProcessor(HDF5Container): """ This class enable to compute and save derived quantities from the raw output @@ -285,6 +302,28 @@ class PostProcessor(HDF5Container): centers = 0.5 * (edges[1:] + edges[:-1]) return (np.stack([values, centers]), {"logbins": logbins}) + def _Brho(self, bins=100): + self.load_cells() + B = _getter_B_int(self.cells) + rho = _getter_rho(self.cells) + B = np.log10(B) + rho = np.log10(rho) + abscisse = np.linspace(np.min(rho), np.max(rho), bins) + weights = mass_func(self.cells) + bin_number = np.zeros(len(B)) + for rho_min in abscisse[:-1]: + bin_number = bin_number + (rho > rho_min).astype(int) + + B_mean = np.zeros(len(abscisse) - 1) + for i in range(len(B_mean)): + B_mean[i] = np.mean(B[bin_number == i]) + values, edges = np.histogram(B, abscisse, weights=weights) + centers = 0.5 * (edges[1:] + edges[:-1]) + print("centers") + if self.pp_params.process.unload_cells: + self.unload_cells() + return {"rho": centers, "B": B_mean} + def _mwa_sigma(self, axes=["x", "y", "z"]): mw_speed = self.save.get_node("/globals/mwa_speed").read() @@ -341,6 +380,15 @@ class PostProcessor(HDF5Container): def _B_v(self, ax_los, z=0.0): return self._vector_v("Br", "unit_mag", ax_los, z) + def _B_int(self, ax_los, z=0.0): + B_op = ScalarOperator( + lambda dset: np.sqrt(np.sum(dset["Br"] ** 2, axis=1)), + self._ro.info["unit_mag"], + ) + dmap_B = (slicing.SliceMap(self._amr, self._cam[ax_los], B_op, z=z)).map.T + dmap_rho = self.save.get_node("/maps/rho_{}".format(ax_los)).read() + return dmap_B + 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 @@ -749,6 +797,14 @@ class PostProcessor(HDF5Container): "/maps", dependencies=["radial_bins", "rr"], ), + "B_int": Rule( + self, + self._B_int, + "Magnetic intensity slice", + "/maps", + dependencies=["rho"], + unit=self.info["unit_mag"], + ), # PDF "rho_pdf": Rule( self, @@ -771,6 +827,20 @@ class PostProcessor(HDF5Container): "/hist", unit=self.info["unit_pressure"], ), + "cos_pdf": Rule( + self, + partial(self._vol_pdf, _getter_abs_cos_vB, logbins=False), + "Global cos-PDF", + "/hist", + unit=cst.none, + ), + "Brho": Rule( + self, + partial(self._Brho), + "Brho average", + "/datasets", + unit=self.info["unit_mag"], + ), # Profiles "axis": Rule( self,