# coding: utf-8 import sys import os import pymses import numpy as np import matplotlib as mpl if os.environ.get("DISPLAY", "") == "": print("No display found. Using non-interactive Agg backend") mpl.use("Agg") import pylab as P import glob as glob try: import cPickle as pickle except: print("cPickle not found, using pickle") import pickle import tables from scipy.stats import linregress 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 # extension for out files out_ext = ".jpeg" P.rcParams["image.cmap"] = "plasma" P.rcParams["savefig.dpi"] = 400 def make_image_disk( path, num, 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"]) rad = 0.5 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, pos_star=np.array([1.0, 1.0, 1.0]), put_title=True, interactive=False, fft=False, ): """ 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_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 """ 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, (-radius + center[1]) * lbox, (radius + center[1]) * lbox, ] time = ro.info["time"] # time in codeunits title = "t=" + str(time)[0:5] + " (code)" # Determining outpout directory if path_out is not None: directory = path_out else: directory = path # Checking for existing files name = directory + "/coldens_z" + "_" + tag + "_" + format(num, "05") + out_ext if len(glob.glob(name)) == 1 and not force: return # Prepare saving data if save_data: maps_disk = {"time": time, "im_extent": im_extent} 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 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] cam = Camera( center=center, line_of_sight_axis=ax_los, region_size=[2.0 * radius, 2.0 * radius], distance=radius, far_cut_depth=radius, up_vector=ax_v, 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: 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.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) 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 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 # 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 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.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.savefig(name_im) P.close() # Toomre parameter if ax_los == "z": def omega_rho_func(dset): pos = dset.get_cell_centers() pos = pos - (pos_star / lbox) xx = pos[:, :, 0] yy = pos[:, :, 1] rc = np.sqrt(xx ** 2 + yy ** 2) # cylindrical radius vx = dset["vel"][:, :, 0] vy = dset["vel"][:, :, 1] omega_rho = 1.0 / rc ** 2 omega_rho = omega_rho * dset["rho"] vyx = vy * xx vxy = vx * yy omega_rho = omega_rho * (vyx - vxy) return omega_rho omega_op = FractionOperator( omega_rho_func, lambda dset: dset["rho"], 1.0 / ro.info["unit_time"] ) cs_op = FractionOperator( lambda dset: dset["P"], lambda dset: dset["rho"], ro.info["unit_velocity"], ) rt_omega = raytracing.RayTracer(amr, ro.info, omega_op) if fft: rt_cs = splatting.SplatterProcessor(amr, ro.info, cs_op, surf_qty=False) else: rt_cs = raytracing.RayTracer(amr, ro.info, cs_op) dmap_omega = rt_omega.process(cam) 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 im = P.imshow( map_Q, extent=im_extent, origin="lower", norm=mpl.colors.LogNorm() ) 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.savefig(name_im) P.close() if cpu: cpu_op = ScalarOperator( lambda dset: dset.icpu * (np.ones(dset["P"].shape)), ro.info["unit_pressure"], ) rt_cpu = raytracing.RayTracer(amr, ro.info, cpu_op) datamap = rt_cpu.process(cam, surf_qty=True) map_cpu = datamap.map.T # if save_data: # 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 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 if interactive: P.figure() else: P.savefig(name_im) P.close() if save_data: name_save = ( directory + "/maps_disk" + "_" + tag + "_" + format(num, "05") + ".save" ) f = open(name_save, "w") pickle.dump(maps_disk, f) f.close() return maps_disk def disk_prop( path_in, num, path_out=None, force=False, nb_bin=20, rad_ext=1.0, mass_star=1.0, pos_star=np.array([1.0, 1.0, 1.0]), binning="log", ): """ Compute properties of a disk in the plane (0,x,y) with a protostar at the center of the box The region of the disk is defined by its scale height Parameters ---------- path_in path of the input data files (output of ramses) num id of the output path_out optional path to the output files force if set, redo ouptut even if already done nb_bin Number of radial bins rad_ext Outer radius of the disk pos_star position of the central protostar mass_star mass of the central protostar """ # Set th output directory if path_out is not None: directory_out = path_out else: directory_out = path_in # Check if the output file exists, and exit if it is the case name_save = directory_out + "/prop_disk_" + str(num).zfill(5) + ".save" if not force and len(glob.glob(name_save)) != 0: return # Compute the bins array if binning == "log": lrad = np.log10(rad_ext) rad = np.logspace(lrad - 2.0, lrad, num=nb_bin) elif binning == "lin": rad = np.linspace(0.0, rad_ext, num=nb_bin) else: raise ValueError( "Incorrect binning specification, binnning should be 'lin' or 'log'" ) # Get Ramses data ro = pymses.RamsesOutput(path_in, num) lbox = ro.info["boxlen"] # boxlen in codeunits (=>pc) time = ro.info["time"] # time in codeunits # Get array of cell positions amr = ro.amr_source(["rho", "vel", "Br", "Bl", "P", "g", "phi"]) cell_source = CellsToPoints(amr) cells = cell_source.flatten() dx = cells.get_sizes() * lbox 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"] # Get radial component of velocity norm_pos = rc norm_pos[rc == 0] = 1.0e-10 # Avoid division per 0 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 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 print("Number of selected cells ", np.sum(mask)) pos_disk = pos[mask] rc_disk = rc[mask] vel_disk = vel[mask] rho_disk = cells["rho"][mask] # density press_disk = cells["P"][mask] # pressure dx_disk = dx[mask] dvol_disk = dx_disk ** 3 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] 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) print("Mass disk", total_mass_disk) print("Mass box", total_mass) # Initialize binned quantities 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) 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) vol_rad = np.zeros(nb_bin - 1) # Volume of a bin surf_rad = np.zeros(nb_bin - 1) # Surface of a bin mass_rad = np.zeros(nb_bin - 1) # Mass of a bin for i in range(nb_bin - 1): 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] rho_rad[i] = np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) / vol_rad[i] temp_rad[i] = np.sum(press_disk[mask_bin] * dvol_disk[mask_bin]) / mass_rad[i] # Surface of a bin : S = dr * 2 * pi * r with # dr = rad[i + 1] - rad[i] and r = (rad[i + 1] + rad[i]) / 2. surf_rad[i] = (rad[i + 1] - rad[i]) * (rad[i + 1] + rad[i]) * np.pi coldens_rad[i] = mass_rad[i] / surf_rad[i] v_az_rad[i] = ( np.sum(v_az_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin]) / mass_rad[i] ) v_rad_rad[i] = ( np.sum(v_rad_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin]) / mass_rad[i] ) height_rad[i] = ( np.sum(height_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin]) / mass_rad[i] ) alpha_rey_rad[i] = (2.0 / 3) * ( ( np.sum( v_az_disk[mask_bin] * v_rad_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin] ) / np.sum(dvol_disk[mask_bin] * press_disk[mask_bin]) - v_az_rad[i] * v_rad_rad[i] * rho_rad[i] / press_rad[i] ) * 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]) / mass_rad[i] ) # Derived quantities cs_rad = np.sqrt(temp_rad) Q_kepl_rad = cs_rad * v_az_rad / (np.pi * G * coldens_rad * rad[0 : nb_bin - 1]) # Means mask_mean = (0.1 < rad[0 : nb_bin - 1]) & (rad[0 : nb_bin - 1] < 0.2) mass_mean = np.sum(mass_rad[mask_mean]) alpha_rey_mean = np.sum(alpha_rey_rad[mask_mean] * mass_rad[mask_mean]) / mass_mean alpha_grav_mean = ( np.sum(alpha_grav_rad[mask_mean] * mass_rad[mask_mean]) / mass_mean ) Q_mean = np.sum(Q_kepl_rad[mask_mean] * mass_rad[mask_mean]) / mass_mean Q_min = np.nanmin(Q_kepl_rad) # store the results prop_disk = { "time": time, "mass_disk": total_mass_disk, "mass_box": total_mass, "rad": rad[:-1], "center": pos_star, "alpha_rey": alpha_rey_rad, "alpha_rey_mean": alpha_rey_mean, "alpha_grav": alpha_grav_rad, "alpha_grav_mean": alpha_grav_mean, "v_rad": v_rad_rad, "v_az": v_az_rad, "v_kepl": v_kepl_rad, "coldens": coldens_rad, "rho": rho_rad, "press": press_rad, "temp": temp_rad, "cs": cs_rad, "Q_kepl": Q_kepl_rad, "Q_mean": Q_mean, "Q_min": Q_min, "height": height_rad, } f = open(name_save, "w") pickle.dump(prop_disk, f) f.close() def plot_disk_prop(path, num, force=False, tag="", interactive=False): """ Plot properties of a disk num id of the ramses output path path to the properties file force if set, redo plots even if already done """ # Load property file name_save = path + "/prop_disk_" + str(num).zfill(5) + ".save" # Check if the properties file exists if len(glob.glob(name_save)) == 0: raise ("no pickle file for disk properties. Run disk_prop() first") f = open(name_save, "r") prop_disk = pickle.load(f) f.close() # Check if the output file exists, and exit if it is the case name_save = path + "/rho_disk_r_" + str(num).zfill(5) + out_ext if not force and len(glob.glob(name_save)) != 0: return 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() 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.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + out_ext) P.close() 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") P.title(title) if interactive: P.figure() else: 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.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 ") P.title(title) if interactive: P.figure() else: 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}$", ) 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.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + out_ext) P.close() # Q P.ylim([0, 10.0]) P.xlim([0, 0.5]) P.yticks(np.arange(0.0, 11, 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 ") P.title(title + ", mass of disk = {} (code)".format(mass)) if interactive: pass else: 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 ") P.title(title + ", mass of box = {} (code)".format(prop_disk["mass_box"])) if interactive: pass else: P.savefig(path + "/H_r_" + str(num).zfill(5) + out_ext) P.close() def disk_pdf( path, num, maps_disk, pos_star=[1.0, 1.0], force=False, interactive=False, nb_bin_hist=50, tag="", rad_min=0.05, rad_max=0.4, put_title=True, ): # Load property file name_prop = path + "/prop_disk_" + str(num).zfill(5) + ".save" # Check if the properties file exists if len(glob.glob(name_prop)) == 0: raise IOError("no pickle file for disk properties. Run disk_prop() first") f = open(name_prop, "r") prop_disk = pickle.load(f) f.close() # Check if the job has already been done if not force: try: slope = prop_disk["fit"]["slope"] print("PDF already computed, slope = {}, exiting ...".format(slope)) return except KeyError: pass # Load maps file print("load maps file") name_maps = path + "/maps_disk" + "_" + tag + "_" + format(num, "05") + ".save" # Check if the maps file exists 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 = prop_disk["time"] 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) + 0.1 ) # radial bins 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)) 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) # Find appropriate bin for each coordinate set bins = np.zeros(rr.shape, dtype=int) for r in rad_bins[1:]: bins = bins + (rr >= r).astype(int) 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_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] # Histogramm : density fluctuation distribution function dcoldens = np.log10(coldens_map.flatten() / coldens_mean) mask_hist = (rr_flat > rad_min) & (rr_flat < rad_max) nb_cells = np.sum(mask_hist) P.grid() P.yscale("log") P.ylim([0.5 / nb_cells, 1.0]) P.xlabel(r"$\log(N / \bar{N})$") P.ylabel("Probability density") if put_title: P.title(title) values, edges, _ = P.hist( dcoldens[mask_hist], nb_bin_hist, range=(-1, 3), weights=np.ones(nb_cells) / nb_cells, ) centers = 0.5 * (edges[1:] + edges[:-1]) # Compute the slope of the right part of the histogramm mask_fit = (centers > 0.25) & (centers < 1.25) & (values > 0) (a, b, rho, _, _) = linregress(centers[mask_fit], np.log10(values[mask_fit])) print("a=%e, b=%e, rho=%e" % (a, b, rho)) fit = { "beta": int(tag.split("_")[1][4:]), "slope": a, "origin": b, "correlation": rho, } P.plot(centers, 10 ** (a * centers + b), "--") f = open(name_prop, "w") prop_disk["fit"] = fit pickle.dump(prop_disk, f) f.close() if interactive: pass else: P.savefig(path + "/dcol_hist_" + tag + "_" + str(num).zfill(5) + out_ext) P.close() # Fluctuation map dcoldens_map = np.reshape(dcoldens, coldens_map.shape) im = P.imshow(dcoldens_map, extent=im_extent, origin="lower", cmap="viridis") if put_title: P.title(title) P.xlabel(r"$x$") P.ylabel(r"$y$") cbar = P.colorbar(im) cbar.set_label(r"$log(N/\bar{N})$ (code)") name = path + "/dcoldensmap_" + tag + "_" + format(num, "05") name_im = name + out_ext if interactive: P.figure() else: P.savefig(name_im) P.close() def compare( path, runs, nums, force=False, interactive=False, Q_in_name=True, pdf=False, gamma=5.0 / 3.0, ): """ Compare time averaged properties of a disk in several simulations nums id or array of ids of the ramses output runs list of runs to consider path path to the properties file force if set, redo plots even if already done interactive interactive mode, to use in a %pylab ipython shell """ if type(nums) == int: nums = [nums] nums_name = "_".join(str(num).zfill(5) for num in nums) # Initialize arrays alpha_rey = np.zeros(len(runs)) alpha_grav = np.zeros(len(runs)) Q = np.zeros(len(runs)) if Q_in_name: Q0 = np.zeros(len(runs)) if pdf: beta = np.zeros(len(runs)) slope = np.zeros(len(runs)) for i, run in enumerate(runs): path_run = path + "/" + run for j, num in enumerate(nums): # Load property file name_save = path_run + "/prop_disk_" + str(num).zfill(5) + ".save" # Check if the properties file exists if len(glob.glob(name_save)) == 0: raise ("no pickle file for disk properties. Run disk_prop() first") f = open(name_save, "r") prop_disk = pickle.load(f) f.close() alpha_rey[i] = alpha_rey[i] + prop_disk["alpha_rey_mean"] / len(nums) alpha_grav[i] = alpha_grav[i] + prop_disk["alpha_grav_mean"] / len(nums) Q[i] = Q[i] + prop_disk["Q_min"] / len(nums) if pdf: fit = prop_disk["fit"] beta[i] = fit["beta"] slope[i] = slope[i] + fit["slope"] / len(nums) if Q_in_name: Q0[i] = float(run.split("_")[2][1:]) # Check if the output file exists, and exit if it is the case name_save = path + "/alphaQ_" + nums_name + out_ext # if (not force and len(glob.glob(name_save)) !=0): # return # alpha = f(Qmin) P.yscale("log") P.ylim([1e-7, 1.0]) P.grid() P.plot(Q, abs(alpha_rey), "o--", label=r"$\bar{\alpha}_{Reynolds}$") P.plot(Q, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$") P.legend() P.ylabel(r"$\bar{\alpha}$") P.xlabel(r"$Q_{min}$") if interactive: P.figure() else: P.savefig(path + "/alphaQ_" + nums_name + out_ext) P.close() if Q_in_name: # alpha = f(Q0) P.yscale("log") P.ylim([1e-7, 1.0]) P.grid() P.plot(Q0, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$") P.plot(Q0, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$") P.legend() P.ylabel(r"$\bar{\alpha}$") P.xlabel(r"$Q_{0}$") if interactive: P.figure() else: P.savefig(path + "/alphaQ0_" + nums_name + out_ext) P.close() if pdf: # slope of the pdf = f(beta) P.grid() P.plot(beta, slope, "o-.") P.legend() P.ylabel(r"$d\log\% / d\logN$") P.xlabel(r"$\beta$") if interactive: P.figure() else: P.savefig(path + "/dcolslopebeta_" + nums_name + out_ext) P.close() # alpha = f(beta) P.yscale("log") P.ylim([1e-7, 1.0]) P.grid() # theoritical alpha (Gammie 2001) alpha_th = (4.0 / 9.0) / (gamma * (gamma - 1.0) * beta) P.plot(beta, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$") P.plot(beta, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$") P.plot(beta, abs(alpha_th), ":", label=r"$\bar{\alpha}_{th}$") P.legend() P.ylabel(r"$\bar{\alpha}$") P.xlabel(r"$\beta_{cool}$") if interactive: P.figure() else: P.savefig(path + "/alphabeta_" + nums_name + out_ext) P.close() def evolution(path, nums, force=False, interactive=False, pdf=False): """ Plot properties over time path path to the properties file nums list of id of the ramses output force if set, redo plots even if already done interactive interactive mode, to use in a %pylab ipython shell """ # Initialize arrays time = np.zeros(len(nums)) alpha_rey = np.zeros(len(nums)) alpha_grav = np.zeros(len(nums)) Qmin = np.zeros(len(nums)) Qmean = np.zeros(len(nums)) mass_disk = np.zeros(len(nums)) mass_box = np.zeros(len(nums)) if pdf: slope = np.zeros(len(nums)) for i, num in enumerate(nums): # Load property file name_prop = path + "/prop_disk_" + str(num).zfill(5) + ".save" # Check if the properties file exists if len(glob.glob(name_prop)) == 0: raise ("no pickle file for disk properties. Run disk_prop() first") f = open(name_prop, "r") prop_disk = pickle.load(f) f.close() time[i] = prop_disk["time"] alpha_rey[i] = prop_disk["alpha_rey_mean"] alpha_grav[i] = prop_disk["alpha_grav_mean"] Qmin[i] = prop_disk["Q_min"] Qmean[i] = prop_disk["Q_mean"] mass_disk[i] = prop_disk["mass_disk"] mass_box[i] = prop_disk["mass_box"] if pdf: slope = prop_disk["fit"]["slope"] # Check if the output file exists, and exit if it is the case name_save = path + "/alpha_time" + out_ext if not force and len(glob.glob(name_save)) != 0: return # Alpha P.yscale("log") P.ylim([1e-7, 1.0]) P.grid() P.plot(time, abs(alpha_rey), "o-.", label=r"$\bar{\alpha}_{Reynolds}$") P.plot(time, abs(alpha_grav), "*--", label=r"$\bar{\alpha}_{grav}$") P.legend() P.ylabel(r"$\bar{\alpha}$") P.xlabel(r"time (code)") if interactive: P.figure() else: P.savefig(path + "/alpha_time" + out_ext) P.close() # Q P.grid() P.plot(time, Qmin, "o-.", label=r"$Q_{min}$") P.plot(time, Qmean, "*--", label=r"$\bar{Q}$") P.legend() P.ylabel(r"$Q$") P.xlabel(r"time (code)") if interactive: P.figure() else: P.savefig(path + "/Q_time" + out_ext) P.close() # M P.grid() P.plot(time, mass_disk, "o-.", label=r"$M_{disk}$") P.plot(time, mass_box, "*--", label=r"$M_{box}$") P.legend() P.ylabel(r"$M / M_{*}$") P.xlabel(r"time (code)") if interactive: P.figure() else: P.savefig(path + "/mass_time" + out_ext) P.close() # PDF if pdf: P.grid() P.plot(time, slope, "o-.") P.legend() P.ylabel(r"$d\log\% / d\logN$") P.xlabel(r"time (code)") if interactive: P.figure() else: P.savefig(path + "/dcolslope_time" + out_ext) P.close()