diff --git a/disk_postprocess.py b/disk_postprocess.py index bdc5bf7..3ef29df 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -24,6 +24,7 @@ from pymses.sources.ramses import output from pymses.analysis import Camera, raytracing, slicing, splatting from pymses.filters import CellsToPoints from pymses.analysis import ScalarOperator, FractionOperator, MaxLevelOperator +from scipy.stats import gaussian_kde # extension for out files out_ext = ".jpeg" @@ -32,84 +33,20 @@ P.rcParams["image.cmap"] = "plasma" P.rcParams["savefig.dpi"] = 400 -def make_image_disk( +def compute_image_data( path, num, - rad=0.5, + radius=0.5, path_out=None, order="<", save_data=True, - force=False, - tag="", - vel_red=20, map_size=512, - put_title=True, - cpu=False, - level=False, - pos_star=np.array([1.0, 1.0, 1.0]), - interactive=False, - fft=False, -): - """ - Make several useful image of an output of a simulation - - Parameters - ---------- - path path of the Ramses output - num Ramses output number - path_out path of the pipeline outputb - order '<' or '>' order used by pymses for reading ramses output - force if set, erase any existing pipeline output files - tag string to add to the output name - vel_red number of point where velocity should be plot in the slices - map_size size of the map - cpuamr plot also levels and cpus at each step - """ - - ro = pymses.RamsesOutput(path, num, order=order) - amr = ro.amr_source(["rho", "vel", "P"]) - center = [0.5, 0.5, 0.5] - - return make_image_aux( - amr, - ro, - center, - rad, - num, - path, - force=force, - path_out=path_out, - save_data=save_data, - map_size=map_size, - vel_red=vel_red, - tag=tag, - cpu=cpu, - level=level, - put_title=put_title, - pos_star=pos_star, - interactive=interactive, - fft=fft, - ) - - -def make_image_aux( - amr, - ro, - center, - radius, - num, - path, - force=False, - path_out=None, - save_data=True, - map_size=512, - vel_red=20, tag="", - cpu=False, - level=False, + axes_los=["x", "y", "z"], + images=["coldens", "rho", "speed", "Q", "T", "level", "cpu"], pos_star=np.array([1.0, 1.0, 1.0]), put_title=True, - interactive=False, + force=False, fft=False, ): """ @@ -117,24 +54,20 @@ def make_image_aux( Parameters ---------- - amr - ro pymses.RamsesOutput object - - center 3D array for coordinates center num output number path_out path of the pipeline output force if set, erase any existing pipeline output files tag string to add to the output name - vel_red number of point where velocity should be plot in the slices map_size size of the map """ + ro = pymses.RamsesOutput(path, num, order=order) + amr = ro.amr_source(["rho", "vel", "P"]) + + center = [0.5, 0.5, 0.5] lbox = ro.info["boxlen"] # boxlen in codeunits G = 1.0 # Gravitational constant - # Plotting parameters - ntick = 6 - title_ax = {"x": "x (code)", "y": "y (code)", "z": "z (code)"} im_extent = [ (-radius + center[0]) * lbox, (radius + center[0]) * lbox, @@ -145,7 +78,7 @@ def make_image_aux( time = ro.info["time"] # time in codeunits title = "t=" + str(time)[0:5] + " (code)" - # Determining outpout directory + # Determining output directory if path_out is not None: directory = path_out else: @@ -157,26 +90,36 @@ def make_image_aux( return # Prepare saving data - if save_data: - maps_disk = {"time": time, "im_extent": im_extent} + if "T" in images and not "rho" in images: + images.append("rho") + if "Q" in images and not "coldens" in images: + images.append("coldens") + maps_disk = { + "time": time, + "im_extent": im_extent, + "center": center, + "radius": radius, + "lbox": lbox, + "images": images, + "axes_los": axes_los, + } + # Prepare raytracing rho_op = ScalarOperator(lambda dset: dset["rho"], ro.info["unit_density"]) - rt = None if fft: rt = splatting.SplatterProcessor(amr, ro.info, rho_op) else: rt = raytracing.RayTracer(amr, ro.info, rho_op) - axes_los = ["x", "y", "z"] # Line of sight axes - axes_h = ["y", "x", "x"] # Horizontal axes - axes_v = ["z", "z", "y"] # Vertical axes + # Prepare axes + ax_nb = {"x": 0, "y": 1, "z": 2} # Associated number of each axes + axes_h = {"x": "y", "y": "x", "z": "x"} # Associated horizontal axe + axes_v = {"x": "z", "y": "z", "z": "y"} # Associated vertical axe - ax_nb = {"x": 0, "y": 1, "z": 2} - - for i, ax_los in enumerate(axes_los): - ax_h = axes_h[i] - ax_v = axes_v[i] + for ax_los in axes_los: + ax_h = axes_h[ax_los] + ax_v = axes_v[ax_los] cam = Camera( center=center, @@ -188,179 +131,40 @@ def make_image_aux( map_max_size=map_size, ) - if interactive: - P.figure() - else: - P.close() - - # Levels - if level: - amr.set_read_levelmax(20) - level_op = MaxLevelOperator() - rt_level = raytracing.RayTracer(amr, ro.info, level_op) - datamap = rt_level.process(cam, surf_qty=True) - map_level = datamap.map.T - # if save_data: - # maps_disk['levels_' + ax_los] = map_level - - levels_ar = np.arange(ro.info["levelmin"], ro.info["levelmax"] + 1) - - # Computing linewidths - lw = np.ones(levels_ar.size) * 2 - lvl_th = 8 # Level threeshold for reducing linewidths - 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=im_extent, - 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) - - if save_data: + if "coldens" in images: + datamap = rt.process(cam, surf_qty=True) + dmap_col = datamap.map.T * lbox maps_disk["coldens_" + ax_los] = dmap_col - im = P.imshow(map_col, extent=im_extent, origin="lower") - - P.locator_params(axis=ax_h, nbins=ntick) - P.locator_params(axis=ax_v, 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"$log(N)$ (code)") - name = directory + "/coldens_" + ax_los + "_" + tag + "_" + format(num, "05") - name_im = name + out_ext - - if interactive: - P.figure() - else: - P.tight_layout(pad=1) - P.savefig(name_im) - P.close() - # Rho slice - datamap_rho = slicing.SliceMap(amr, cam, rho_op, z=0.0) - dmap_rho = (datamap_rho.map).T - map_rho = np.log10(dmap_rho) + if "rho" in images: + datamap_rho = slicing.SliceMap(amr, cam, rho_op, z=0.0) + dmap_rho = (datamap_rho.map).T + maps_disk["rho_" + ax_los] = dmap_rho - vh_op = ScalarOperator( - lambda dset: dset["vel"][:, ax_nb[ax_h]], ro.info["unit_velocity"] - ) - dmap_vh = slicing.SliceMap(amr, cam, vh_op, z=0.0) - map_vh_red = dmap_vh.map[ - ::vel_red, ::vel_red - ] # take only a subset of velocities - map_vh_red = map_vh_red.T + if "speed" in images: + vh_op = ScalarOperator( + lambda dset: dset["vel"][:, ax_nb[ax_h]], ro.info["unit_velocity"] + ) + dmap_vh = slicing.SliceMap(amr, cam, vh_op, z=0.0) + vv_op = ScalarOperator( + lambda dset: dset["vel"][:, ax_nb[ax_v]], ro.info["unit_velocity"] + ) + dmap_vv = slicing.SliceMap(amr, cam, vv_op, z=0.0) - vv_op = ScalarOperator( - lambda dset: dset["vel"][:, ax_nb[ax_v]], ro.info["unit_velocity"] - ) - dmap_vv = slicing.SliceMap(amr, cam, vv_op, z=0.0) - map_vv_red = dmap_vv.map[::vel_red, ::vel_red] - map_vv_red = map_vv_red.T + maps_disk["v" + ax_h + "_" + ax_los] = (dmap_vh.map).T + maps_disk["v" + ax_v + "_" + ax_los] = (dmap_vv.map).T - # if save_data: - # maps_disk['rho_' + ax_los] = dmap_rho - # maps_disk['v' + ax_h + '_' + ax_los] = (dmap_vh.map).T - # maps_disk['v' + ax_v + '_' + ax_los] = (dmap_vv.map).T + if "T" in images: + P_op = ScalarOperator(lambda dset: dset["P"], ro.info["unit_pressure"]) + dmap_P = (slicing.SliceMap(amr, cam, P_op, z=0.0)).map.T - im = P.imshow(map_rho, extent=im_extent, origin="lower") - P.locator_params(axis=ax_h, nbins=ntick) - P.locator_params(axis=ax_v, nbins=ntick) - nh = map_vh_red.shape[0] - nv = map_vv_red.shape[1] - vec_h = ( - np.arange(nh) * 2.0 / nh * radius - radius + center[0] + radius / nh - ) * lbox - vec_v = ( - np.arange(nv) * 2.0 / nv * radius - radius + center[1] + radius / nv - ) * lbox - hh, vv = np.meshgrid(vec_h, vec_v) - max_v = np.max(np.sqrt(map_vh_red ** 2 + map_vv_red ** 2)) - - Q = P.quiver(hh, vv, map_vh_red, map_vv_red, units="width") - P.quiverkey( - Q, - 0.7, - 0.95, - max_v, - r"$" + str(max_v)[0:4] + "$ (code)", - labelpos="E", - coordinates="figure", - ) - 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"$log(n)$ (code)") - - name = directory + "/rho_" + ax_los + "_" + tag + "_" + format(num, "05") - name_im = name + out_ext - if interactive: - P.figure() - else: - P.tight_layout(pad=1) - P.savefig(name_im) - P.close() - - P_op = ScalarOperator(lambda dset: dset["P"], ro.info["unit_pressure"]) - dmap_P = (slicing.SliceMap(amr, cam, P_op, z=0.0)).map.T - - dmap_T = dmap_P / dmap_rho - map_T = np.log10(dmap_T) - - # if save_data: - # maps_disk['T_' + ax_los] = dmap_T - - im = P.imshow(map_T, extent=im_extent, 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"$log(T) \, (K)$") - - name = directory + "/T_" + ax_los + "_" + tag + "_" + format(num, "05") - name_im = name + out_ext - if interactive: - P.figure() - else: - P.tight_layout(pad=1) - P.savefig(name_im) - P.close() + dmap_T = dmap_P / dmap_rho + maps_disk["T_" + ax_los] = dmap_T # Toomre parameter - - if ax_los == "z": + if "Q" in images and ax_los == "z": def omega_rho_func(dset): pos = dset.get_cell_centers() @@ -397,39 +201,19 @@ def make_image_aux( dmap_cs = rt_cs.process(cam) map_Q = (lbox * dmap_cs.map.T) * dmap_omega.map.T / (np.pi * G * dmap_col) - # if save_data: - # maps_disk['Q_' + ax_los] = map_Q + maps_disk["Q_" + ax_los] = map_Q - im = P.imshow( - map_Q, - extent=im_extent, - origin="lower", - cmap="RdBu", - norm=mpl.colors.LogNorm(), - vmin=0.01, - vmax=100.0, - ) + # Levels + if "levels" in images: + amr.set_read_levelmax(20) + level_op = MaxLevelOperator() + rt_level = raytracing.RayTracer(amr, ro.info, level_op) + datamap = rt_level.process(cam, surf_qty=True) + map_level = datamap.map.T + maps_disk["levels_" + ax_los] = map_level - 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"$Q$") - - name = directory + "/Q_" + ax_los + "_" + tag + "_" + format(num, "05") - name_im = name + out_ext - if interactive: - P.figure() - else: - P.tight_layout(pad=1) - P.savefig(name_im) - P.close() - - if cpu: + # Cpus + if "cpu" in images: cpu_op = ScalarOperator( lambda dset: dset.icpu * (np.ones(dset["P"].shape)), ro.info["unit_pressure"], @@ -438,35 +222,228 @@ def make_image_aux( datamap = rt_cpu.process(cam, surf_qty=True) map_cpu = datamap.map.T - # if save_data: - # maps_disk['cpu_' + ax_los] = map_cpu + maps_disk["cpu_" + ax_los] = map_cpu - im = P.imshow(map_cpu, extent=im_extent, origin="lower") - P.locator_params(axis="x", nbins=ntick) - P.locator_params(axis="y", nbins=ntick) + if save_data: + f = open(name_save, "w") + pickle.dump(maps_disk, f) + f.close() + return maps_disk + + +def plot_maps( + path, + num, + force=False, + vel_red=40, + tag="", + images=None, + axes_los=None, + maps_disk=None, + interactive=False, + put_title=True, +): + """ + Make several useful image of an output of a simulation, auxillary function + + Parameters + ---------- + amr + ro pymses.RamsesOutput object + + center 3D array for coordinates center + num output number + path path of the pipeline output + force if set, erase any existing pipeline output files + tag string to add to the output name + vel_red Take point each vel_red for velocities + map_size size of the map + """ + + path_out = path + + # Load maps file + print("load maps file") + name_maps = path + "/maps_disk" + "_" + tag + "_" + format(num, "05") + ".save" + + if maps_disk is None: + if len(glob.glob(name_maps)) == 0: + raise IOError("no pickle file for disk maps. Run make_image_disk() first") + f = open(name_maps, "r") + maps_disk = pickle.load(f) + f.close() + print("maps file loaded") + + time = maps_disk["time"] + im_extent = maps_disk["im_extent"] + center = maps_disk["center"] + radius = maps_disk["radius"] + lbox = maps_disk["lbox"] + + # Plot parameters + title = "t=" + str(time)[0:5] + " (code)" + ntick = 6 + title_ax = {"x": r"$x$ (code)", "y": r"$y$ (code)", "z": r"$z$ (code)"} + + # Determine output directory + if path_out is not None: + directory = path_out + else: + directory = path + + # Determine what to plot + if images == None: + images = maps_disk["images"] + + if axes_los == None: + axes_los = maps_disk["axes_los"] + + # Prepare axes + axes_h = {"x": "y", "y": "x", "z": "x"} # Associated horizontal axe + axes_v = {"x": "z", "y": "z", "z": "y"} # Associated vertical axe + + P.close() + + for ax_los in axes_los: + ax_h = axes_h[ax_los] + ax_v = axes_v[ax_los] + + for image in images: + + if ( + (image == "Q" and not ax_los == "z") + or image == "levels" + or image == "speed" + ): + continue + + map_disk = maps_disk[image + "_" + ax_los] + + if image == "Q": + im = P.imshow( + map_Q, + extent=im_extent, + origin="lower", + cmap="RdBu", + norm=mpl.colors.LogNorm(), + vmin=0.01, + vmax=100.0, + ) + else: + im = P.imshow( + map_disk, + extent=im_extent, + origin="lower", + norm=mpl.colors.LogNorm(), + ) + + P.locator_params(axis=ax_h, nbins=ntick) + P.locator_params(axis=ax_v, 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"cpu") - name = directory + "/cpu_" + ax_los + "_" + tag + "_" + format(num, "05") - name_im = name + out_ext + cbar = P.colorbar(im) + + if image == "coldens": + cbar.set_label(r"$log(N)$ (code)") + + if "levels" in images: + map_level = maps_disk["levels_" + ax_los] + # Computing linewidths + lw = np.ones(levels_ar.size) * 2 + lvl_th = 8 # Level threeshold for reducing linewidths + 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=im_extent, + 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", + ) + elif image == "rho": + cbar.set_label(r"$log(n)$ (code)") + + if "speed" in images: + dmap_vh = maps_disk["v" + ax_h + "_" + ax_los] + dmap_vv = maps_disk["v" + ax_v + "_" + ax_los] + + map_vh_red = dmap_vh[ + ::vel_red, ::vel_red + ] # take only a subset of velocities + map_vv_red = dmap_vv[::vel_red, ::vel_red] + nh = map_vh_red.shape[0] + nv = map_vv_red.shape[1] + vec_h = ( + np.arange(nh) * 2.0 / nh * radius + - radius + + center[0] + + radius / nh + ) * lbox + vec_v = ( + np.arange(nv) * 2.0 / nv * radius + - radius + + center[1] + + radius / nv + ) * lbox + hh, vv = np.meshgrid(vec_h, vec_v) + max_v = np.max(np.sqrt(map_vh_red ** 2 + map_vv_red ** 2)) + + Q = P.quiver(hh, vv, map_vh_red, map_vv_red, units="width") + P.quiverkey( + Q, + 0.7, + 0.95, + max_v, + r"$" + str(max_v)[0:4] + "$ (code)", + labelpos="E", + coordinates="figure", + ) + + elif image == "T": + cbar.set_label(r"$log(T) \, (K)$") + elif image == "Q": + cbar.set_label(r"$Q$") + else: + cbar.set_label(image) + + name = ( + directory + + "/" + + image + + "_" + + ax_los + + "_" + + tag + + "_" + + format(num, "05") + + out_ext + ) if interactive: P.figure() else: P.tight_layout(pad=1) - P.savefig(name_im) + P.savefig(name) P.close() - if save_data: - f = open(name_save, "w") - pickle.dump(maps_disk, f) - f.close() return maps_disk @@ -502,7 +479,7 @@ def disk_prop( mass_star mass of the central protostar """ - # Set th output directory + # Set the output directory if path_out is not None: directory_out = path_out else: @@ -559,8 +536,10 @@ def disk_prop( 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_dens = cells["rho"] > 0.01 * np.mean(cells["rho"]) # condition on density + mask_vel = abs(v_rad / v_az) < 1.0 # condition on speed + + mask = mask_pos # & mask_dens & mask_vel print("Number of selected cells ", np.sum(mask)) pos_disk = pos[mask] @@ -573,7 +552,7 @@ def disk_prop( v_rad_disk = v_rad[mask] v_az_disk = v_az[mask] v_kepl = np.sqrt(mass_star * G / rc_disk) - height_disk = height[mask] + # height_disk = height[mask] g_rad_disk = g_rad[mask] g_az_disk = g_az[mask] @@ -605,7 +584,6 @@ def disk_prop( mask_bin = (rc_disk > rad[i]) & (rc_disk < rad[i + 1]) # print("Bin #{} : {} cells between {} and {}".format(i, np.sum(mask_bin), rad[i], rad[i + 1])) - vol_rad[i] = np.sum(dvol_disk[mask_bin]) mass_rad[i] = np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) press_rad[i] = np.sum(press_disk[mask_bin] * dvol_disk[mask_bin]) / vol_rad[i] @@ -627,10 +605,12 @@ def disk_prop( / mass_rad[i] ) - height_rad[i] = ( - np.sum(height_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin]) - / mass_rad[i] - ) + try: + height_rad[i] = ( + np.max(pos_disk[:, 2][mask_bin]) - np.min(pos_disk[:, 2][mask_bin]) + ) / 2.0 + except ValueError: + height_rad[i] = 0 alpha_rey_rad[i] = (2.0 / 3) * ( ( @@ -707,7 +687,15 @@ def disk_prop( f.close() -def plot_disk_prop(path, num, force=False, tag="", interactive=False, put_title=False): +def plot_disk_prop( + path, + num, + plots=["rho", "T", "V", "coldens", "Q", "alpha", "H"], + force=False, + tag="", + interactive=False, + put_title=False, +): """ Plot properties of a disk @@ -733,154 +721,99 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False, put_title= time = prop_disk["time"] mass = prop_disk["mass_disk"] - title = "t=" + str(time)[0:5] + " (code)" rad = prop_disk["rad"] - if interactive: - P.figure() - else: - P.close() + for plot in plots: + title = "t=" + str(time)[0:5] + " (code)" + P.grid() + P.xlabel("disk radius") + if plot == "rho": + P.xscale("log") + P.yscale("log") + P.plot(rad, prop_disk["rho"], color="k", linewidth=2) + P.ylabel(r"$n \, (code)$") + elif plot == "T": + P.xscale("log") + P.yscale("log") + P.plot(rad, prop_disk["temp"], color="k", linewidth=2) + P.ylabel(r"$T \, (K)$") + elif plot == "V": + P.xscale("log") + P.yscale("symlog", linthreshy=0.01) + P.plot(rad, prop_disk["v_rad"], color="k", linewidth=2) + P.plot(rad, prop_disk["v_kepl"], color="b", linewidth=2) + P.plot(rad, abs(prop_disk["v_az"]), color="r", linewidth=2) + P.plot(rad, prop_disk["cs"], color="c", linewidth=2) + P.legend( + (r"$v_r$", r"$v_{kepl}$", r"$v_\phi$", r"$c_s$"), loc="upper right" + ) + P.ylabel(r"$V \, (km s^{-1})$") + elif plot == "coldens": + P.xscale("log") + P.yscale("log") + P.plot(rad, prop_disk["coldens"], color="k", linewidth=2) + P.ylabel(r"$N\, (cm^{-2})$") + elif plot == "alpha": + alpha_rey_mean, alpha_grav_mean = ( + prop_disk["alpha_rey_mean"], + prop_disk["alpha_grav_mean"], + ) + P.xlim([1e-2, 0.25]) + P.yscale("log") + P.ylim([1e-7, 1.0]) + P.plot( + rad, + abs(prop_disk["alpha_rey"]), + "b", + linewidth=2, + label=r"$\alpha_{Reynolds}$", + ) + P.plot(rad, abs(alpha_rey_mean * np.ones(len(rad))), "b:", linewidth=1) + P.plot( + rad, + abs(prop_disk["alpha_grav"]), + "r", + linewidth=2, + label=r"$\alpha_{grav}$", + ) + P.plot(rad, abs(alpha_grav_mean * np.ones(len(rad))), "r:", linewidth=1) + P.plot( + rad, + abs(prop_disk["alpha_rey"]) + abs(prop_disk["alpha_grav"]), + "g--", + linewidth=2, + label=r"$\alpha_{tot}$", + ) + P.legend() + P.ylabel(r"$\alpha$") + title = ( + title + + r", $\bar{\alpha}_{Reynolds} = %.1e, \bar{\alpha}_{grav} = %.1e$" + % (alpha_rey_mean, alpha_grav_mean) + ) + elif plot == "Q": + P.ylim([0, 3.0]) + P.xlim([0, 0.25]) + P.yticks(np.arange(0.0, 3, 1.0)) + P.plot(rad, abs(prop_disk["Q_kepl"]), color="b", linewidth=2) + # P.plot(rad, abs(prop_disk['Q_mean']) * np.ones(len(rad)), 'b:', linewidth=1) + P.ylabel(r"$Q$") + P.xlabel("disk radius ") + title = title + ", mass of disk = {} (code)".format(mass) + elif plot == "H": + P.plot(rad, abs(prop_disk["height"] / rad), color="b", linewidth=2) + P.ylabel(r"H ratio") + title = title + ", mass of box = {} (code)".format(prop_disk["mass_box"]) - P.xscale("log") - P.yscale("log") - P.grid() - P.plot(rad, prop_disk["rho"], color="k", linewidth=2) - P.ylabel(r"$n \, (code)$") - P.xlabel("disk radius") - P.title(title) - if interactive: - P.figure() - else: - P.tight_layout(pad=1) - P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + out_ext) - P.close() + if put_title: + P.title(title) - P.xscale("log") - P.yscale("log") - P.grid() - P.plot(rad, prop_disk["temp"], color="k", linewidth=2) - P.ylabel(r"$T \, (K)$") - P.xlabel("disk radius") - if put_title: - P.title(title) - if interactive: - P.figure() - else: - P.tight_layout(pad=1) - P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + out_ext) - P.close() - - P.xscale("log") - P.yscale("symlog", linthreshy=0.01) - - P.plot(rad, prop_disk["v_rad"], color="k", linewidth=2) - P.plot(rad, prop_disk["v_kepl"], color="b", linewidth=2) - P.plot(rad, abs(prop_disk["v_az"]), color="r", linewidth=2) - P.plot(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})$") - P.xlabel("disk radius") - if interactive: - P.figure() - else: - P.tight_layout(pad=1) - P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + out_ext) - P.close() - - P.xscale("log") - P.yscale("log") - P.grid() - P.plot(rad, prop_disk["coldens"], color="k", linewidth=2) - P.ylabel(r"$N\, (cm^{-2})$") - P.xlabel("disk radius ") - if put_title: - P.title(title) - if interactive: - P.figure() - else: - P.tight_layout(pad=1) - P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + out_ext) - P.close() - - # Alpha - alpha_rey_mean, alpha_grav_mean = ( - prop_disk["alpha_rey_mean"], - prop_disk["alpha_grav_mean"], - ) - - # P.xscale('log') - P.xlim([1e-2, 0.25]) - P.yscale("log") - P.ylim([1e-7, 1.0]) - P.grid() - P.plot( - rad, abs(prop_disk["alpha_rey"]), "b", linewidth=2, label=r"$\alpha_{Reynolds}$" - ) - P.plot(rad, abs(alpha_rey_mean * np.ones(len(rad))), "b:", linewidth=1) - P.plot( - rad, abs(prop_disk["alpha_grav"]), "r", linewidth=2, label=r"$\alpha_{grav}$" - ) - P.plot(rad, abs(alpha_grav_mean * np.ones(len(rad))), "r:", linewidth=1) - P.plot( - rad, - abs(prop_disk["alpha_rey"]) + abs(prop_disk["alpha_grav"]), - "g--", - linewidth=2, - label=r"$\alpha_{tot}$", - ) - if put_title: - P.title( - title - + r", $\bar{\alpha}_{Reynolds} = %.1e, \bar{\alpha}_{grav} = %.1e$" - % (alpha_rey_mean, alpha_grav_mean) - ) - - P.legend() - P.ylabel(r"$\alpha$") - P.xlabel("disk radius ") - - if interactive: - P.figure() - else: - P.tight_layout(pad=1) - P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + out_ext) - P.close() - - # Q - P.ylim([0, 3.0]) - P.xlim([0, 0.25]) - P.yticks(np.arange(0.0, 3, 1.0)) - P.grid() - P.plot(rad, abs(prop_disk["Q_kepl"]), color="b", linewidth=2) - # P.plot(rad, abs(prop_disk['Q_mean']) * np.ones(len(rad)), 'b:', linewidth=1) - P.ylabel(r"$Q$") - P.xlabel("disk radius ") - if put_title: - P.title(title + ", mass of disk = {} (code)".format(mass)) - - if interactive: - pass - else: - P.tight_layout(pad=1) - P.savefig(path + "/Q_r_" + str(num).zfill(5) + out_ext) - P.close() - - # height ratio - P.grid() - P.plot(rad, abs(prop_disk["height"] / rad), color="b", linewidth=2) - P.ylabel(r"H ratio") - P.xlabel("disk radius ") - if put_title: - P.title(title + ", mass of box = {} (code)".format(prop_disk["mass_box"])) - - if interactive: - pass - else: - P.tight_layout(pad=1) - P.savefig(path + "/H_r_" + str(num).zfill(5) + out_ext) - P.close() + if interactive: + P.figure() + else: + P.tight_layout(pad=1) + P.savefig(path + "/" + plot + "_disk_r_" + str(num).zfill(5) + out_ext) + P.close() def disk_pdf( @@ -893,8 +826,9 @@ def disk_pdf( nb_bin_hist=50, tag="", rad_min=0.05, - rad_max=0.4, + rad_max=0.3, put_title=True, + do_speed=True, ): # Load property file @@ -930,10 +864,9 @@ def disk_pdf( print("maps file loaded") time = prop_disk["time"] + im_extent = maps_disk["im_extent"] title = tag.split("_")[1] + " t=" + str(time)[0:5] + " (code)" - coldens_map = maps_disk["coldens_z"] - im_extent = maps_disk["im_extent"] # radius of the corner of the box plus a margin rad_box = ( np.sqrt((im_extent[1] - pos_star[0]) ** 2 + (im_extent[3] - pos_star[1]) ** 2) @@ -943,15 +876,28 @@ def disk_pdf( rad_bins = prop_disk["rad"] rad_bins = 0.5 * (rad_bins[0:-1] + rad_bins[1:]) rad_bins = np.concatenate(([0.0], rad_bins, [rad_box])) - coldens_rad = prop_disk["coldens"] - # Add value for borders - coldens_rad = np.concatenate(([coldens_rad[0]], coldens_rad)) + + coldens_map = maps_disk["coldens_z"] x = np.linspace(im_extent[0], im_extent[1], coldens_map.shape[0]) y = np.linspace(im_extent[2], im_extent[3], coldens_map.shape[1]) xx, yy = np.meshgrid(x, y) rr = np.sqrt((xx - pos_star[0]) ** 2 + (yy - pos_star[1]) ** 2) + rho_map = maps_disk["rho_z"] + cs_map = np.sqrt(maps_disk["T_z"]) + + vx_map = maps_disk["vx_z"] + vy_map = maps_disk["vy_z"] + v_map = np.sqrt(vx_map ** 2 + vy_map ** 2) + v_kepl = np.sqrt(1.0 / rr) + xx_star = xx - pos_star[0] + yy_star = yy - pos_star[1] + vrad_map = (vx_map * xx_star + vy_map * yy_star) / rr + vaz_map = (vy_map * xx_star - vx_map * yy_star) / rr + + maps = [coldens_map, rho_map, cs_map, v_map, vaz_map] + # Find appropriate bin for each coordinate set bins = np.zeros(rr.shape, dtype=int) for r in rad_bins[1:]: @@ -959,19 +905,35 @@ def disk_pdf( bins_flat = bins.flatten() rr_flat = rr.flatten() - # Retrieve the mean of the bins - # use linear interpolation to improve accuracy - coldens_mean = (rad_bins[bins_flat + 1] - rr_flat) * coldens_rad[bins_flat] - coldens_mean = ( - coldens_mean + (rr_flat - rad_bins[bins_flat]) * coldens_rad[bins_flat + 1] - ) - coldens_mean = coldens_mean / (rad_bins[bins_flat + 1] - rad_bins[bins_flat]) + coldens_rad = np.zeros(len(rad_bins) - 1) + rho_rad = np.zeros(len(rad_bins) - 1) + cs_rad = np.zeros(len(rad_bins) - 1) + v_rad = np.zeros(len(rad_bins) - 1) + vaz_rad = np.zeros(len(rad_bins) - 1) + rads = [coldens_rad, rho_rad, cs_rad, v_rad, vaz_rad] + means = [None] * len(rads) + + for i, rad_arr in enumerate(rads): + map_arr = maps[i] + + for j in range(len(rad_bins) - 1): + rad_arr[j] = np.mean(map_arr[bins == j]) + + # Add value for borders + rad_arr = np.concatenate(([rad_arr[0]], rad_arr)) + + # Retrieve the mean of the bins + # use linear interpolation to improve accuracy + means[i] = (rad_bins[bins_flat + 1] - rr_flat) * rad_arr[bins_flat] + means[i] = means[i] + (rr_flat - rad_bins[bins_flat]) * rad_arr[bins_flat + 1] + means[i] = means[i] / (rad_bins[bins_flat + 1] - rad_bins[bins_flat]) + + coldens_mean = means[0] coldens_mean_map = np.reshape(coldens_mean, coldens_map.shape) - # Kill fluctuations for r < rad_min and r > rad_max - coldens_map[rr < rad_min] = coldens_mean_map[rr < rad_min] - coldens_map[rr > rad_max] = coldens_mean_map[rr > rad_max] + mask_map = (rr > rad_min) & (rr < rad_max) + coldens_map[np.logical_not(mask_map)] = np.nan # Histogramm : density fluctuation distribution function dcoldens = np.log10(coldens_map.flatten() / coldens_mean) @@ -992,6 +954,8 @@ def disk_pdf( ) centers = 0.5 * (edges[1:] + edges[:-1]) + # Variance + var = np.var(dcoldens[mask_hist]) # Compute the slope of the right part of the histogramm mask_fit = (centers > 0.25) & (centers < 1.25) & (values > 0) if np.sum(mask_fit > 0): @@ -999,13 +963,18 @@ def disk_pdf( centers[mask_fit], np.log10(values[mask_fit]) ) P.plot(centers, 10 ** (a * centers + b), "--", linewidth=2) - print("a=%e, b=%e, rho=%e" % (a, b, rho)) + print("a=%e, b=%e, rho=%e, var=%e" % (a, b, rho, var)) + try: + beta = int(tag.split("_")[1][4:]) + except ValueError: + beta = 0 fit = { - "beta": int(tag.split("_")[1][4:]), + "beta": beta, "slope": a, "origin": b, "correlation": rho, "stderr": stderr, + "var": var, } f = open(name_prop, "w") prop_disk["fit"] = fit @@ -1013,13 +982,65 @@ def disk_pdf( f.close() if interactive: - pass + P.figure() else: P.tight_layout(pad=1) P.savefig(path + "/dcol_hist_" + tag + "_" + str(num).zfill(5) + out_ext) P.close() - # Fluctuation map + # Fluctuations plots + rho_mean = means[1] + cs_mean = means[2] + v_mean = means[3] + vaz_mean = means[4] + drho = abs(rho_map.flatten() / rho_mean) + dcs = abs(cs_map.flatten() / cs_mean) + dvaz = abs(vaz_map.flatten() / vaz_mean) + dv = abs(v_map.flatten() / v_mean) + dmach = abs((v_map.flatten() - v_mean) / cs_map.flatten()) + dmach_mean = np.mean(dmach[mask_hist]) + f = open(name_prop, "w") + prop_disk["dmach_mean"] = dmach_mean + pickle.dump(prop_disk, f) + f.close() + + # P.scatter(drho[mask_hist], dcs[mask_hist], s=0.1) + P.hist2d( + np.log10(drho[mask_hist]), + np.log10(dcs[mask_hist]), + (1000, 1000), + norm=mpl.colors.LogNorm(), + ) + P.xlabel(r"$\log(\rho / \bar{\rho})$") + P.ylabel(r"$\log(c_s / \bar{c_s})$") + P.ylim([-0.6, 0.6]) + + P.colorbar() + if interactive: + P.figure() + else: + P.tight_layout(pad=1) + P.savefig(path + "/dcs_" + tag + "_" + str(num).zfill(5) + out_ext) + P.close() + + P.hist2d( + np.log10(drho[mask_hist]), + dvaz[mask_hist], + (1000, 1000), + norm=mpl.colors.LogNorm(), + ) + P.xlabel(r"$\log(\rho / \bar{\rho})$") + P.ylabel(r"$\log(v_\varphi / \bar{v_\varphi})$") + + P.colorbar() + if interactive: + P.figure() + else: + P.tight_layout(pad=1) + P.savefig(path + "/dvaz_" + tag + "_" + str(num).zfill(5) + out_ext) + P.close() + + # Fluctuations maps dcoldens_map = np.reshape(dcoldens, coldens_map.shape) im = P.imshow(dcoldens_map, extent=im_extent, origin="lower", cmap="viridis") @@ -1041,6 +1062,144 @@ def disk_pdf( P.savefig(name_im) P.close() + dcs_map = np.reshape(dcs, coldens_map.shape) + dcs_map[np.logical_not(mask_map)] = np.nan + + im = P.imshow( + np.log10(dcs_map), + extent=im_extent, + origin="lower", + cmap="RdBu_r", + vmin=-0.6, + vmax=0.6, + ) + + if put_title: + P.title(title) + P.xlabel(r"$x$") + P.ylabel(r"$y$") + + cbar = P.colorbar(im) + cbar.set_label(r"$log(c_s/\bar{c_s})$ (code)") + name = path + "/dcsmap_" + tag + "_" + format(num, "05") + name_im = name + out_ext + + if interactive: + P.figure() + else: + P.tight_layout(pad=1) + P.savefig(name_im) + P.close() + + drho_map = np.reshape(drho, coldens_map.shape) + drho_map[np.logical_not(mask_map)] = np.nan + + im = P.imshow( + np.log10(drho_map), + extent=im_extent, + origin="lower", + cmap="RdBu_r", + vmin=-2.0, + vmax=2.0, + ) + + if put_title: + P.title(title) + P.xlabel(r"$x$") + P.ylabel(r"$y$") + + cbar = P.colorbar(im) + cbar.set_label(r"$log(\rho/\bar{\rho})$ (code)") + name = path + "/drhomap_" + tag + "_" + format(num, "05") + name_im = name + out_ext + + if interactive: + P.figure() + else: + P.tight_layout(pad=1) + P.savefig(name_im) + P.close() + + if do_speed: + dvel_kepl = abs(v_map - v_kepl) / v_kepl + dvaz_kepl = abs(vaz_map - v_kepl) / v_kepl + dvrad_vaz = abs(vrad_map / vaz_map) + dvaz_map = np.reshape(dvaz, coldens_map.shape) + dvaz_map[np.logical_not(mask_map)] = np.nan + + im = P.imshow( + dvaz_kepl, + extent=im_extent, + origin="lower", + cmap="RdBu_r", + norm=mpl.colors.LogNorm(), + vmax=1.0, + vmin=0.01, + ) + cbar = P.colorbar(im) + cbar.set_label(r"$|v_\varphi - v_{kepl}|/v_{kepl}$") + P.xlabel(r"$x$") + P.ylabel(r"$y$") + + if put_title: + P.title(title) + + name = path + "/dvaz_kepl_" + tag + "_" + format(num, "05") + name_im = name + out_ext + + if interactive: + P.figure() + else: + P.tight_layout(pad=1) + P.savefig(name_im) + P.close() + + im = P.imshow( + vrad_map, + extent=im_extent, + origin="lower", + cmap="RdBu_r", + norm=mpl.colors.SymLogNorm(linthresh=0.03), + vmax=8.0, + vmin=-8, + ) + cbar = P.colorbar(im) + cbar.set_label(r"$v_r$") + P.xlabel(r"$x$") + P.ylabel(r"$y$") + + if put_title: + P.title(title) + + name = path + "/vradmap_" + tag + "_" + format(num, "05") + name_im = name + out_ext + + if interactive: + P.figure() + else: + P.tight_layout(pad=1) + P.savefig(name_im) + P.close() + + im = P.imshow(dvaz_map, extent=im_extent, origin="lower", cmap="RdBu_r") + cbar = P.colorbar(im) + cbar.set_label(r"$v_\varphi / \bar{v_\varphi}$") + P.xlabel(r"$x$") + P.ylabel(r"$y$") + + if put_title: + P.title(title) + + name = path + "/dvazmap_" + tag + "_" + format(num, "05") + name_im = name + out_ext + + if interactive: + P.figure() + else: + P.tight_layout(pad=1) + P.savefig(name_im) + P.close() + def compare( path, @@ -1078,12 +1237,18 @@ def compare( beta = np.zeros(len(runs)) slope = np.zeros(len(runs)) slope_std = np.zeros(len(runs)) + var = np.zeros(len(runs)) + var_std = np.zeros(len(runs)) + dmach = np.zeros(len(runs)) + dmach_std = np.zeros(len(runs)) for i, run in enumerate(runs): path_run = path + "/" + run nb_outputs = 0 slope_run = [] + var_run = [] + dmach_run = [] for num in nums: try: @@ -1108,6 +1273,8 @@ def compare( fit = prop_disk["fit"] beta[i] = fit["beta"] slope_run.append(fit["slope"]) + var_run.append(fit["var"]) + dmach_run.append(prop_disk["dmach_mean"]) nb_outputs = nb_outputs + 1 print(run, num, nb_outputs) @@ -1122,11 +1289,17 @@ def compare( if pdf: slope[i] = np.mean(slope_run) slope_std[i] = np.std(slope_run) + var[i] = np.mean(var_run) + var_std[i] = np.std(var_run) + dmach[i] = np.mean(dmach_run) + dmach_std[i] = np.std(dmach_run) else: for array in [alpha_rey, alpha_grav, Q]: array[i] = np.nan if pdf: slope[i] = np.nan + var[i] = np.nan + dmach[i] = np.nan if Q_in_name: Q0[i] = float(run.split("_")[2][1:]) @@ -1195,6 +1368,44 @@ def compare( P.savefig(path_out + "/dcolslopebeta_" + nums_name + out_ext) P.close() + # var of the pdf = f(beta) + P.grid() + P.errorbar(beta, var, yerr=var_std, fmt="o") + + P.legend() + P.ylabel(r"$Var(\log(N / \bar(N))$") + P.xlabel(r"$\beta$") + + (a, b, rho, _, stderr) = linregress(beta, var) + P.plot(beta, a * beta + b, "--", linewidth=1.5) + print("a=%e, b=%e, rho^2=%e" % (a, b, rho ** 2)) + + if interactive: + P.figure() + else: + P.tight_layout(pad=1) + P.savefig(path_out + "/varbeta_" + nums_name + out_ext) + P.close() + + # dmach = f(var) + P.grid() + P.errorbar(var, dmach, xerr=var_std, yerr=dmach_std, fmt="o") + + P.legend() + P.ylabel(r"$<(v - \bar{v}) / c_s>$") + P.xlabel(r"$Var(\log(N / \bar(N))$") + + (a, b, rho, _, stderr) = linregress(var, dmach) + P.plot(var, a * var + b, "--", linewidth=1.5) + print("a=%e, b=%e, rho^2=%e" % (a, b, rho ** 2)) + + if interactive: + P.figure() + else: + P.tight_layout(pad=1) + P.savefig(path_out + "/dmachvar_" + nums_name + out_ext) + P.close() + # alpha = f(beta) P.yscale("log") P.ylim([1e-7, 1.0]) diff --git a/pipeline_disk.py b/pipeline_disk.py index 95caa80..f26bf1b 100644 --- a/pipeline_disk.py +++ b/pipeline_disk.py @@ -50,7 +50,8 @@ parser.add_argument( parser.add_argument( "-pd", "--plot_disk", help="plot specific disk radial analysis", action="store_true" ) -parser.add_argument("-m", "--maps", help="do generic maps", action="store_true") +parser.add_argument("-m", "--maps", help="compute generic maps", action="store_true") +parser.add_argument("-pm", "--plot_maps", help="plot generic maps", action="store_true") parser.add_argument( "-c", "--compare", help="compare different runs", action="store_true" ) @@ -67,8 +68,15 @@ parser.add_argument( parser.add_argument( "--fft", help="use quick and dirty fft rendering", action="store_true" ) -parser.add_argument("--level", help="plot levels", action="store_true") -parser.add_argument("--cpu", help="plot cpu", action="store_true") +parser.add_argument( + "--images", + nargs="*", + default=["coldens", "rho", "speed", "Q", "T"], + choices=["coldens", "rho", "speed", "Q", "T", "level", "cpu"], +) +parser.add_argument( + "--axes", nargs="*", default=["x", "y", "z"], choices=["x", "y", "z"] +) parser.add_argument("--zoom", help="zoom", type=float, default=2.0) parser.add_argument( "-ms", @@ -78,36 +86,6 @@ parser.add_argument( default=1024, ) -parser.add_argument( - "-d", "--disk", help="do specific disk radial analysis", action="store_true" -) -parser.add_argument("-m", "--maps", help="do generic maps", action="store_true") -parser.add_argument( - "-c", "--compare", help="compare different runs", action="store_true" -) -parser.add_argument( - "-ev", - "--evolution", - help="plot evolution of quantities over time", - action="store_true", -) -parser.add_argument( - "--pdf", help="plot pdf of fluctuations of column density", action="store_true" -) - -parser.add_argument( - "--fft", help="use quick and dirty fft rendering", action="store_true" -) -parser.add_argument("--level", help="plot levels", action="store_true") -parser.add_argument("--cpu", help="plot cpu", action="store_true") -parser.add_argument( - "-ms", - "--mapsize", - help="size of the maps created in he map mode (in pixel)", - type=int, - default=1024, -) - parser.add_argument( "--nb_bin", help="Number of bins for azimuthal averages", type=int, default=50 @@ -200,22 +178,34 @@ for run in runs: maps_disk = None if args.maps: print("[{}, {}] computing maps".format(run, i)) - maps_disk = dp.make_image_disk( + maps_disk = dp.compute_image_data( path_in, i, - rad=rad, + radius=rad, path_out=path_out, tag=run, map_size=args.mapsize, force=args.force_redo, - level=args.level, - cpu=args.cpu, + axes_los=args.axes, + images=args.images, pos_star=np.array([args.x, args.y, args.z]), fft=args.fft, + ) + print("[{}, {}] maps computed".format(run, i)) + if args.plot_maps: + print("[{}, {}] plotting maps".format(run, i)) + maps_disk = dp.plot_maps( + path_out, + i, + maps_disk=maps_disk, + axes_los=args.axes, + images=args.images, + tag=run, + force=args.force_redo, interactive=args.interactive, put_title=(not args.beamer), ) - print("[{}, {}] maps computed".format(run, i)) + print("[{}, {}] maps plotted".format(run, i)) if args.disk: print("[{}, {}] computing disk props".format(run, i)) dp.disk_prop(