diff --git a/disk_postprocess.py b/disk_postprocess.py index 0ba4cef..ef99e73 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -9,7 +9,6 @@ import pickle as pickle import module_extract as me from pymses.filters import CellsToPoints - from pymses.sources.ramses import output from pymses.analysis import Camera, raytracing, slicing @@ -343,7 +342,7 @@ def disk_prop( path_out=None, force=False, nb_bin=20, - rad_ext=100.0, + rad_ext=1.0, mass_star=1.0, pos_star=np.array([1.0, 1.0, 1.0]), ): @@ -387,42 +386,21 @@ def disk_prop( ro = pymses.RamsesOutput(path_in, num) lbox = ro.info["boxlen"] # boxlen in codeunits (=>pc) - ( - AU, - pc, - Ms, - Myr, - scale_n, - scale_d, - scale_t, - scale_l, - scale_v, - scale_T2, - scale_ener, - scale_mag, - microG, - km_s, - Cwnm, - scale_mass, - unit_col, - lbox_pc, - ) = me.normalisation(ro) - - time = ro.info["time"] * scale_t / Myr + time = ro.info["time"] # * scale_t / Myr # Get array of cell positions amr = ro.amr_source(["rho", "vel", "Br", "Bl", "P"]) cell_source = CellsToPoints(amr) cells = cell_source.flatten() dx = cells.get_sizes() - pos = cells.points + pos = cells.points * lbox # Get positions in the frame of the protostar pos = pos - pos_star # Get cylindrical radius rc = np.sqrt(pos[:, 0] ** 2 + pos[:, 1] ** 2) # Get velocities - vel = cells["vel"] + vel = cells["vel"] * lbox # Get radial component of velocity norm_pos = rc norm_pos[rc == 0] = 1.0e-10 # Avoid division per 0 @@ -432,12 +410,12 @@ def disk_prop( # Select cells that are actually in the disk, ie within the scale height # TODO Check units - G = 6.8e-8 - cs = np.sqrt(cells["P"] / cells["rho"]) * scale_v # sound velocity + G = 1.0 # G=6.8e-8 + cs = np.sqrt(cells["P"] / cells["rho"]) # sound velocity height = cs * np.sqrt(rc ** 3 / (G * mass_star)) mask_pos = np.abs(pos[:, 2]) < height # condition on position mask_dens = cells["rho"] > 1.0e6 # condition on density - mask = mask_pos | mask_dens + mask = mask_pos # & mask_dens print("Number of selected cells ", np.sum(mask)) pos_disk = pos[mask] @@ -449,33 +427,25 @@ def disk_prop( dvol_disk = dx_disk ** 3 v_rad_disk = v_rad[mask] v_az_disk = v_az[mask] - - # TODO Check what do that does - nzoom = 9 - eps = 0.5 ** nzoom - # map_coldens , map_w13, xedges, yedges = me.make_hierarch_map(pos_disk_x,pos_disk_y,pos_disk_z,dx_disk,rho_disk,rho_disk,eps,center=[0.,0.,0.],make_image=do_plot,path_out=directory_out,tag='xy_'+ str(num).zfill(5)) - - # map_coldens , map_w13, xedges, yedges = me.make_hierarch_map(pos_disk_z,pos_disk_x,pos_disk_y,dx_disk,rho_disk,rho_disk,eps,center=[0.,0.,0.],make_image=do_plot,path_out=directory_out,tag='xz_'+ str(num).zfill(5)) + v_kepl = np.sqrt(mass_star * G / rc_disk) # Initialize binned quantities - norm_rad = lbox * scale_l / AU # radius in AU - rdisk_AU = rc_disk * norm_rad - cs_rad = np.zeros(nb_bin - 1) temp_rad = np.zeros(nb_bin - 1) press_rad = np.zeros(nb_bin - 1) rho_rad = np.zeros(nb_bin - 1) coldens_rad = np.zeros(nb_bin - 1) v_az_rad = np.zeros(nb_bin - 1) + v_kepl_rad = np.zeros(nb_bin - 1) v_rad_rad = np.zeros(nb_bin - 1) alpha_rey_rad = np.zeros(nb_bin - 1) for i in range(nb_bin - 1): - mask_bin = (rdisk_AU > rad[i]) & (rdisk_AU < rad[i + 1]) + mask_bin = (rc_disk > rad[i]) & (rc_disk < rad[i + 1]) print( - "Bin {} cells between {} and {} AU".format( - np.sum(mask_bin), rad[i], rad[i + 1] + "Bin #{} : {} cells between {} and {}".format( + i, np.sum(mask_bin), rad[i], rad[i + 1] ) ) @@ -492,8 +462,8 @@ def disk_prop( # 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]) - * (lbox * pc) ** 3 - / ((rad[i + 1] - rad[i]) * (rad[i + 1] + rad[i]) * np.pi * AU ** 2) + * (lbox) ** 3 + / ((rad[i + 1] - rad[i]) * (rad[i + 1] + rad[i]) * np.pi) ) v_az_rad[i] = np.sum( @@ -519,13 +489,18 @@ def disk_prop( / abs(v_az_rad[i]) ) - # 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_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]) - v_az_rad = v_az_rad * scale_v / km_s - v_rad_rad = v_rad_rad * scale_v / km_s + # 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 prop_disk = { "time": time, @@ -534,6 +509,7 @@ def disk_prop( "alpha_rey": alpha_rey_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, @@ -577,7 +553,7 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): linewidth=2, ) P.ylabel(r"$\log(n) \, (cm^{-3})$") - P.xlabel("disk radius (AU)") + P.xlabel("disk radius") if pdf: P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + ".pdf") @@ -591,7 +567,7 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): linewidth=2, ) P.ylabel(r"$\log(T) \, (K)$") - P.xlabel("disk radius (AU)") + P.xlabel("disk radius") if pdf: P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + ".pdf") @@ -603,27 +579,19 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): P.yscale("symlog", linthreshy=0.01) P.plot((prop_disk["rad_AU"]), ((prop_disk["v_rad"])), color="k", linewidth=2) + P.plot((prop_disk["rad_AU"]), ((prop_disk["v_kepl"])), color="b", linewidth=2) P.plot((prop_disk["rad_AU"]), (abs(prop_disk["v_az"])), color="r", linewidth=2) P.plot((prop_disk["rad_AU"]), ((prop_disk["cs"])), color="c", linewidth=2) - P.legend((r"$v_r$", r"$v_\phi$", r"$c_s$"), loc="upper right") + P.legend((r"$v_r$", r"$v_{kepl}$", r"$v_\phi$", r"$c_s$"), loc="upper right") P.ylabel(r"$V \, (km s^{-1})$") - P.xlabel("disk radius (AU)") + P.xlabel("disk radius") if pdf: P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + ".pdf") P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + ".jpeg") - P.legend((r"$v_r$", r"$v_\phi$", r"$c_s$"), loc="upper right") - - P.ylabel(r"$V \, (km s^{-1})$") - P.xlabel("disc radius (AU)") - - if pdf: - P.savefig(path + "V_disk_r_" + str(num).zfill(5) + ".pdf") - P.savefig(path + "V_disk_r_" + str(num).zfill(5) + ".jpeg") - P.clf() P.plot( np.log10(prop_disk["rad_AU"]), @@ -632,7 +600,7 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): linewidth=2, ) P.ylabel(r"$\log(N) \, (cm^{-2})$") - P.xlabel("disk radius (AU)") + P.xlabel("disk radius ") if pdf: P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + ".pdf") @@ -642,12 +610,12 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): P.xscale("log") P.yscale("symlog", linthreshy=0.001) - P.plot(prop_disk["rad_AU"], prop_disk["alpha_rey"], color="b", linewidth=2) + P.plot(prop_disk["rad_AU"], abs(prop_disk["alpha_rey"]), color="b", linewidth=2) - P.plot(prop_disc["rad_AU"], prop_disc["alpha_rey"], color="b", linewidth=2) + # P.legend(r'$\alpha_{Rey}$', loc='upper right') - P.ylabel(r"$\alpha}$") - P.xlabel("disk radius (AU)") + P.ylabel(r"$\alpha$") + P.xlabel("disk radius ") if pdf: P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + ".pdf") diff --git a/pipeline_disk.py b/pipeline_disk.py index eedc592..058087b 100644 --- a/pipeline_disk.py +++ b/pipeline_disk.py @@ -21,6 +21,10 @@ parser.add_argument( "-l", "--last_output", help="id of last output", type=int, default=100 ) parser.add_argument("-s", "--step", help="step between two output", type=int, default=1) +parser.add_argument( + "-d", "--disk", help="do specific disk radial analysis", action="store_true" +) + args = parser.parse_args() @@ -54,5 +58,9 @@ for run in runs: mag_im=False, AU_units=False, ) - dp.disk_prop(path_in, i, path_out=path_out, rad_ext=50000) - dp.plot_disk_prop(path_out, i, tag=run + "_") + # me.look(path_in, i) + if args.disk: + dp.disk_prop( + path_in, i, path_out=path_out, rad_ext=1, nb_bin=50, force=True + ) + dp.plot_disk_prop(path_out, i, tag=run + "_")