Overplot levels, add pdf of density fluctuations

disk_postprocess.py

@@ -37,6 +37,8 @@ def make_image_disk(
     map_size=512,
     put_title=True,
     cpuamr=False,
+    cpu=False,
+    level=False,
     pos_star=np.array([1.0, 1.0, 1.0]),
     interactive=False,
     fft=False,
@@ -75,6 +77,8 @@ def make_image_disk(
         vel_red=vel_red,
         tag=tag,
         cpuamr=cpuamr,
+        cpu=cpu,
+        level=level,
         put_title=put_title,
         pos_star=pos_star,
         interactive=interactive,
@@ -95,6 +99,8 @@ def make_image_aux(
     vel_red=20,
     tag="",
     cpuamr=False,
+    cpu=False,
+    level=False,
     pos_star=np.array([1.0, 1.0, 1.0]),
     put_title=True,
     interactive=False,
@@ -118,6 +124,9 @@ def make_image_aux(
     cpuamr    plot also levels and cpus at each step
     """
 
+    cpu = cpu or cpuamr
+    level = level or cpuamr
+
     lbox = ro.info["boxlen"]  # boxlen in codeunits
     lbox_units = lbox
 
@@ -168,17 +177,55 @@ def make_image_aux(
             map_max_size=map_size,
         )
 
-        datamap = rt.process(cam, surf_qty=True)
-
-        # Column density
-        dmap_col = datamap.map.T * lbox
-        map_col = np.log10(dmap_col)
-
         if interactive:
             P.figure()
         else:
             P.close()
 
+        # Levels
+        if level:
+            level_op = MaxLevelOperator()
+            amr.set_read_levelmax(20)
+            rt_level = raytracing.RayTracer(amr, ro.info, level_op)
+            datamap = rt_level.process(cam, surf_qty=True)
+            map_level = datamap.map.T
+            levels_ar = np.arange(ro.info["levelmin"], ro.info["levelmax"] + 1)
+            # Computing linewidths
+            lw = np.ones(levels_ar.size) * 2
+            lvl_th = 8
+            lw[levels_ar >= lvl_th] = lw[levels_ar >= lvl_th] ** (
+                lvl_th - levels_ar[levels_ar >= lvl_th]
+            )
+            lw[levels_ar < lvl_th] = 1.0
+
+            cont = P.contour(
+                map_level,
+                extent=[
+                    (-radius + center[0]) * lbox_units,
+                    (radius + center[0]) * lbox_units,
+                    (-radius + center[1]) * lbox_units,
+                    (radius + center[1]) * lbox_units,
+                ],
+                origin="lower",
+                colors="k",
+                linewidths=lw,
+                levels=levels_ar,
+            )
+            cont.levels = cont.levels + 1
+
+            P.clabel(
+                cont,
+                levels_ar[levels_ar < lvl_th + 2][1::2],
+                inline=1,
+                fontsize=8.0,
+                fmt="%1d",
+            )
+
+        # Column density
+        datamap = rt.process(cam, surf_qty=True)
+        dmap_col = datamap.map.T * lbox
+        map_col = np.log10(dmap_col)
+
         im = P.imshow(
             map_col,
             extent=[
@@ -385,43 +432,7 @@ def make_image_aux(
                 P.savefig(name_im)
                 P.close()
 
-        if cpuamr:
+        if cpu:
-            level_op = MaxLevelOperator()
-            amr.set_read_levelmax(20)
-            rt_level = raytracing.RayTracer(amr, ro.info, level_op)
-            datamap = rt_level.process(cam, surf_qty=True)
-            map_level = datamap.map.T
-
-            im = P.imshow(
-                map_level,
-                extent=[
-                    (-radius + center[0]) * lbox_units,
-                    (radius + center[0]) * lbox_units,
-                    (-radius + center[1]) * lbox_units,
-                    (radius + center[1]) * lbox_units,
-                ],
-                origin="lower",
-            )
-
-            P.locator_params(axis="x", nbins=ntick)
-            P.locator_params(axis="y", nbins=ntick)
-
-            if put_title:
-                P.title(title)
-            P.xlabel(title_ax[ax_h])
-            P.ylabel(title_ax[ax_v])
-            cbar = P.colorbar(im)
-            cbar.set_label(r"level")
-
-            name = directory + "/level_" + ax_los + "_" + tag + "_" + format(num, "05")
-            name_im = name + out_ext
-
-            if interactive:
-                P.figure()
-            else:
-                P.savefig(name_im)
-                P.close()
-
             cpu_op = ScalarOperator(
                 lambda dset: dset.icpu * (np.ones(dset["P"].shape)),
                 ro.info["unit_pressure"],
@@ -503,6 +514,8 @@ def disk_prop(
     if not force and len(glob.glob(name_save)) != 0:
         return
 
+    nb_bin_hist = nb_bin
+
     # Compute the bins array
     lrad = np.log10(rad_ext)
     rad = np.logspace(lrad - 2.0, lrad, num=nb_bin)
@@ -514,7 +527,7 @@ def disk_prop(
     time = ro.info["time"]  # time in codeunits
 
     # Get array of cell positions
-    amr = ro.amr_source(["rho", "vel", "Br", "Bl", "P"])
+    amr = ro.amr_source(["rho", "vel", "Br", "Bl", "P", "g", "phi"])
     cell_source = CellsToPoints(amr)
     cells = cell_source.flatten()
     dx = cells.get_sizes() * lbox
@@ -532,6 +545,10 @@ def disk_prop(
     v_rad = (pos[:, 0] * vel[:, 0] + pos[:, 1] * vel[:, 1]) / norm_pos
     # Get azimuthal component of velocity
     v_az = (pos[:, 0] * vel[:, 1] - pos[:, 1] * vel[:, 0]) / norm_pos
+    # Gravitational field
+    g = cells["g"]
+    g_rad = (pos[:, 0] * g[:, 0] + pos[:, 1] * g[:, 1]) / norm_pos
+    g_az = (pos[:, 0] * g[:, 1] - pos[:, 1] * g[:, 0]) / norm_pos
 
     # Select cells that are actually in the disk, ie within the scale height
     G = 1.0
@@ -553,6 +570,8 @@ def disk_prop(
     v_az_disk = v_az[mask]
     v_kepl = np.sqrt(mass_star * G / rc_disk)
     height_disk = height[mask]
+    g_rad_disk = g_rad[mask]
+    g_az_disk = g_az[mask]
 
     total_mass_disk = np.sum(rho_disk * dvol_disk)
     total_mass = np.sum(cells["rho"] * dx ** 3)
@@ -570,9 +589,15 @@ def disk_prop(
     v_kepl_rad = np.zeros(nb_bin - 1)
     v_rad_rad = np.zeros(nb_bin - 1)
     alpha_rey_rad = np.zeros(nb_bin - 1)
+    alpha_rey_rad_bis = np.zeros(nb_bin - 1)
+    alpha_grav_rad = np.zeros(nb_bin - 1)
     Q_kepl_rad = np.zeros(nb_bin - 1)
     height_rad = np.zeros(nb_bin - 1)
 
+    # Density fluctuations
+    hist_drho = np.zeros(nb_bin_hist)
+    hist_edges = np.zeros(nb_bin_hist + 1)
+
     for i in range(nb_bin - 1):
         mask_bin = (rc_disk > rad[i]) & (rc_disk < rad[i + 1])
 
@@ -588,7 +613,6 @@ def disk_prop(
             rho_disk[mask_bin] * dvol_disk[mask_bin]
         )
 
-        # TODO verifier unites
         # Surface of a bin :  S = dr * 2 * pi * r with
         #  dr = rad[i + 1] - rad[i] and r = (rad[i + 1] + rad[i]) / 2.
         coldens_rad[i] = np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) / (
@@ -621,19 +645,39 @@ def disk_prop(
             / abs(v_az_rad[i])
         )
 
+        # alpha_rey_rad_bis[i] = (2./3) * (np.sum((v_az_disk[mask_bin] - v_az_rad[i])
+        #                                         * (v_rad_disk[mask_bin] - v_rad_rad[i])
+        #                                         * rho_disk[mask_bin] * dvol_disk[mask_bin])
+        #                                  / np.sum(dvol_disk[mask_bin] * press_disk[mask_bin])
+        #                                  * v_az_rad[i] / abs(v_az_rad[i]))
+
+        alpha_grav_rad[i] = (2.0 / 3) * (
+            np.sum(
+                g_az_disk[mask_bin]
+                * g_rad_disk[mask_bin]
+                * rho_disk[mask_bin]
+                * dvol_disk[mask_bin]
+            )
+            / (4 * np.pi * G)
+            / np.sum(dvol_disk[mask_bin] * press_disk[mask_bin])
+            * coldens_rad[i]
+        )
+
         v_kepl_rad[i] = np.sum(
             v_kepl[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin]
         ) / np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin])
 
-    # Convert to good units (TODO check)
+        # Histogramm : density fluctuaction distribution function
-    cs_rad = np.sqrt(temp_rad)  # *scale_v / km_s
+        drho = np.log(rho_disk[mask_bin] / rho_rad[i])
-    temp_rad = temp_rad  # * scale_T2
+        hist, hist_edges = P.histogram(
-    press_rad = press_rad  # * scale_v**2 * scale_d
+            drho,
-
+            bins=nb_bin_hist,
-    v_az_rad = v_az_rad  # * scale_v / km_s
+            weights=dvol_disk[mask_bin] * 2.0 ** (3 * ro.info["levelmax"]),
-    v_rad_rad = v_rad_rad  # * scale_v / km_s
+        )
-    v_kepl_rad = v_kepl_rad
+        hist_drho = hist_drho + hist
 
+    print(hist_drho, hist_edges)
+    cs_rad = np.sqrt(temp_rad)
     Q_kepl_rad = cs_rad * v_az_rad / (np.pi * G * coldens_rad * rad[0 : nb_bin - 1])
 
     prop_disk = {
@@ -643,12 +687,16 @@ def disk_prop(
         "rad": rad[0 : nb_bin - 1],
         "center": pos_star,
         "alpha_rey": alpha_rey_rad,
+        # 'alpha_rey_bis':alpha_rey_rad_bis,
+        "alpha_grav": alpha_grav_rad,
         "v_rad": v_rad_rad,
         "v_az": v_az_rad,
         "v_kepl": v_kepl_rad,
         "coldens": coldens_rad,
         "rho": rho_rad,
         "press": press_rad,
+        "hist_drho": hist_drho,
+        "hist_edges": hist_edges,
         "temp": temp_rad,
         "cs": cs_rad,
         "Q_kepl": Q_kepl_rad,
@@ -727,7 +775,7 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
     P.plot((prop_disk["rad"]), ((prop_disk["v_kepl"])), color="b", linewidth=2)
     P.plot((prop_disk["rad"]), (abs(prop_disk["v_az"])), color="r", linewidth=2)
     P.plot((prop_disk["rad"]), ((prop_disk["cs"])), color="c", linewidth=2)
-
+    P.grid()
     P.legend((r"$v_r$", r"$v_{kepl}$", r"$v_\phi$", r"$c_s$"), loc="upper right")
 
     P.ylabel(r"$V \, (km s^{-1})$")
@@ -752,14 +800,33 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
         P.close()
 
     # Alpha
-    P.xscale("log")
+    # P.xscale('log')
+    P.xlim([1e-2, 0.25])
     P.yscale("log")
-    P.ylim([1e-5, 1.0])
+    P.ylim([1e-7, 1.0])
-
+    P.grid()
-    P.plot(prop_disk["rad"], abs(prop_disk["alpha_rey"]), color="b", linewidth=2)
+    P.plot(
-
+        prop_disk["rad"],
-    P.plot(prop_disk["rad"], abs(prop_disk["alpha_rey"]), color="b", linewidth=2)
+        abs(prop_disk["alpha_rey"]),
+        linewidth=2,
+        label=r"$\alpha_{Reynolds}$",
+    )
+    # P.plot(prop_disk['rad'],abs(prop_disk['alpha_rey_bis']), '--', linewidth=1,label=r"$\alpha_R 2$")
+    P.plot(
+        prop_disk["rad"],
+        abs(prop_disk["alpha_grav"]),
+        linewidth=2,
+        label=r"$\alpha_{grav}$",
+    )
+    P.plot(
+        prop_disk["rad"],
+        abs(prop_disk["alpha_rey"]) + abs(prop_disk["alpha_grav"]),
+        "--",
+        linewidth=2,
+        label=r"$\alpha_{tot}$",
+    )
 
+    P.legend()
     P.ylabel(r"$\alpha$")
     P.xlabel("disk radius ")
     P.title(title)
@@ -786,7 +853,7 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
         P.savefig(path + "/Q_r_" + str(num).zfill(5) + out_ext)
         P.close()
 
-    # height ration
+    # height ratio
     P.grid()
     P.plot(
         prop_disk["rad"],
@@ -803,3 +870,20 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
     else:
         P.savefig(path + "/H_r_" + str(num).zfill(5) + out_ext)
         P.close()
+
+    # Density fluctuation histogram
+    P.grid()
+    P.xlabel(r"$\log(\frac{\rho}{\bar{\rho}})$")
+    P.ylabel(r"# of cells")
+    P.title(title)
+    hist = prop_disk["hist_drho"]
+    egdes = prop_disk["hist_edges"]
+    widths = egdes[1:] - egdes[:-1]
+    centers = egdes[:-1] + widths / 2.0
+    P.bar(centers, hist, width=widths)
+
+    if interactive:
+        pass
+    else:
+        P.savefig(path + "/drho_hist_" + str(num).zfill(5) + out_ext)
+        P.close()