From 9f070425f4907df85e1e52b75a439b626080bc63 Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Wed, 8 May 2019 12:15:57 +0200 Subject: [PATCH] Overplot levels, add pdf of density fluctuations --- disk_postprocess.py | 206 +++++++++++++++++++++++++++++++------------- 1 file changed, 145 insertions(+), 61 deletions(-) diff --git a/disk_postprocess.py b/disk_postprocess.py index 9a0592c..96fe87f 100644 --- a/disk_postprocess.py +++ b/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: - 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() - + if cpu: 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) - cs_rad = np.sqrt(temp_rad) # *scale_v / km_s - temp_rad = temp_rad # * scale_T2 - press_rad = press_rad # * scale_v**2 * scale_d - - v_az_rad = v_az_rad # * scale_v / km_s - v_rad_rad = v_rad_rad # * scale_v / km_s - v_kepl_rad = v_kepl_rad + # Histogramm : density fluctuaction distribution function + drho = np.log(rho_disk[mask_bin] / rho_rad[i]) + hist, hist_edges = P.histogram( + drho, + bins=nb_bin_hist, + weights=dvol_disk[mask_bin] * 2.0 ** (3 * ro.info["levelmax"]), + ) + 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.plot(prop_disk["rad"], abs(prop_disk["alpha_rey"]), color="b", linewidth=2) - - P.plot(prop_disk["rad"], abs(prop_disk["alpha_rey"]), color="b", linewidth=2) + P.ylim([1e-7, 1.0]) + P.grid() + P.plot( + prop_disk["rad"], + 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()