From db17181fea3b094ea90865261406300bd4bd722b Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Mon, 13 May 2019 16:03:46 +0200 Subject: [PATCH] Add pdf for coldens --- disk_postprocess.py | 386 +++++++++++++++++++++++++++++++++++++------- pipeline_disk.py | 54 ++++++- 2 files changed, 373 insertions(+), 67 deletions(-) diff --git a/disk_postprocess.py b/disk_postprocess.py index e01da7b..2dd2045 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -1,9 +1,9 @@ # coding: utf-8 import sys -import numpy as np import os -import pymses +import pymses +import numpy as np import matplotlib as mpl if os.environ.get("DISPLAY", "") == "": @@ -12,7 +12,13 @@ if os.environ.get("DISPLAY", "") == "": import pylab as P import glob as glob -import pickle as pickle + +try: + import cPickle as pickle +except: + print("cPickle not found, using pickle") + import pickle +import tables from pymses.sources.ramses import output from pymses.analysis import Camera, raytracing, slicing, splatting @@ -31,12 +37,12 @@ def make_image_disk( num, path_out=None, order="<", + save_data=True, force=False, tag="", vel_red=20, map_size=512, put_title=True, - cpuamr=False, cpu=False, level=False, pos_star=np.array([1.0, 1.0, 1.0]), @@ -51,7 +57,7 @@ def make_image_disk( path path of the Ramses output num Ramses output number path_out path of the pipeline outputb - order '<' or '>' TODO + 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 @@ -64,7 +70,7 @@ def make_image_disk( rad = 0.5 center = [0.5, 0.5, 0.5] - make_image_aux( + return make_image_aux( amr, ro, center, @@ -73,10 +79,10 @@ def make_image_disk( path, force=force, path_out=path_out, + save_data=save_data, map_size=map_size, vel_red=vel_red, tag=tag, - cpuamr=cpuamr, cpu=cpu, level=level, put_title=put_title, @@ -95,10 +101,10 @@ def make_image_aux( path, force=False, path_out=None, + save_data=True, map_size=512, vel_red=20, tag="", - cpuamr=False, cpu=False, level=False, pos_star=np.array([1.0, 1.0, 1.0]), @@ -121,39 +127,39 @@ def make_image_aux( 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 """ - cpu = cpu or cpuamr - level = level or cpuamr - lbox = ro.info["boxlen"] # boxlen in codeunits - lbox_units = lbox - 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_units, - (radius + center[0]) * lbox_units, - (-radius + center[1]) * lbox_units, - (radius + center[1]) * lbox_units, + (-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 @@ -189,15 +195,19 @@ def make_image_aux( # Levels if level: - level_op = MaxLevelOperator() 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 + 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] ) @@ -226,6 +236,9 @@ def make_image_aux( 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) @@ -249,9 +262,9 @@ def make_image_aux( P.close() # Rho slice - dmap_rho = slicing.SliceMap(amr, cam, rho_op, z=0.0) - map_rho = np.log10(dmap_rho.map) - map_rho = map_rho.T + 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"] @@ -260,7 +273,6 @@ def make_image_aux( 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( @@ -270,6 +282,11 @@ def make_image_aux( 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) @@ -277,12 +294,13 @@ def make_image_aux( nv = map_vv_red.shape[1] vec_h = ( np.arange(nh) * 2.0 / nh * radius - radius + center[0] + radius / nh - ) * lbox_units + ) * lbox vec_v = ( np.arange(nv) * 2.0 / nv * radius - radius + center[1] + radius / nv - ) * lbox_units + ) * 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, @@ -293,7 +311,6 @@ def make_image_aux( labelpos="E", coordinates="figure", ) - if put_title: P.title(title) P.xlabel(title_ax[ax_h]) @@ -310,11 +327,14 @@ def make_image_aux( P.close() P_op = ScalarOperator(lambda dset: dset["P"], ro.info["unit_pressure"]) - dmap_P = slicing.SliceMap(amr, cam, P_op, z=0.0) + dmap_P = (slicing.SliceMap(amr, cam, P_op, z=0.0)).map.T - dmap_T = dmap_P.map.T / dmap_rho.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) @@ -372,8 +392,10 @@ def make_image_aux( dmap_omega = rt_omega.process(cam) dmap_cs = rt_cs.process(cam) - dmap_Q = (lbox * dmap_cs.map.T) * dmap_omega.map.T / (np.pi * G * dmap_col) - map_Q = dmap_Q + 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() @@ -406,6 +428,9 @@ 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 + im = P.imshow(map_cpu, extent=im_extent, origin="lower") P.locator_params(axis="x", nbins=ntick) P.locator_params(axis="y", nbins=ntick) @@ -427,6 +452,15 @@ def make_image_aux( 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, @@ -437,6 +471,7 @@ def disk_prop( 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) @@ -470,11 +505,16 @@ 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) + 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) @@ -553,10 +593,6 @@ def disk_prop( surf_rad = np.zeros(nb_bin - 1) # Surface of a bin mass_rad = np.zeros(nb_bin - 1) # Mass of a bin - # 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]) @@ -620,20 +656,12 @@ def disk_prop( / mass_rad[i] ) - # 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] - ) - hist_drho = hist_drho + hist - # 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) - print(rad[0 : nb_bin - 1][mask_mean]) 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 = ( @@ -641,7 +669,6 @@ def disk_prop( ) Q_mean = np.sum(Q_kepl_rad[mask_mean] * mass_rad[mask_mean]) / mass_mean Q_min = np.nanmin(Q_kepl_rad) - print("alphas, Q ", alpha_rey_mean, alpha_grav_mean, Q_mean) # store the results prop_disk = { @@ -660,8 +687,6 @@ def disk_prop( "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, @@ -837,21 +862,105 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False): P.savefig(path + "/H_r_" + str(num).zfill(5) + out_ext) P.close() - # Density fluctuation histogram + +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.1, +): + + # 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() + + # 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 = "t=" + str(time)[0:5] + " (code)" + + coldens_map = maps_disk["coldens_z"] + im_extent = maps_disk["im_extent"] + rad_bins = prop_disk["rad"] + coldens_rad = prop_disk["coldens"] + + 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) + + coldens_mean = coldens_rad[bins.flatten()] + coldens_mean_map = np.reshape(coldens_mean, coldens_map.shape) + + # Kill fluctuations for r < rad_min + coldens_map[rr < rad_min] = coldens_mean_map[rr < rad_min] + + # Histogramm : density fluctuation distribution function + dcoldens = np.log10(coldens_map.flatten() / coldens_mean) P.grid() - P.xlabel(r"$\log(\frac{\rho}{\bar{\rho}})$") - P.ylabel(r"fraction of total volume") + P.yscale("log") + P.xlabel(r"$\log(N / \bar{N})$") + P.ylabel("Probability density") 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) + P.hist(dcoldens, nb_bins_hist, range=(-1, 2)) if interactive: pass else: - P.savefig(path + "/drho_hist_" + str(num).zfill(5) + out_ext) + P.savefig(path + "/dcol_hist_" + str(num).zfill(5) + out_ext) + P.close() + + # Fluctuation map + dcoldens_map = np.reshape(dcoldens, coldens_map.shape) + + # cont = P.contour(coldens_mean_map, + # extent=im_extent, + # origin='lower', + # colors='k', + # linewidths=0.1) + + im = P.imshow(dcoldens_map, extent=im_extent, origin="lower", cmap="viridis") + + 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_" + format(num, "05") + name_im = name + out_ext + + if interactive: + P.figure() + else: + P.savefig(name_im) P.close() @@ -1026,3 +1135,162 @@ def evolution(path, nums, force=False, interactive=False): else: P.savefig(path + "/mass_time" + out_ext) P.close() + + +def make_clump_hop( + path_in, + num, + name, + thres_dens, + thres_level, + pos_zoom, + size_zoom, + path_out=None, + path_hop="", + force=False, + gcomp=True, +): + """ + This routine use the HOP algorithm to extract clumps defined from their peaks + as an output it provides a list of cell position ordered by the group to which they belong + + Parameters + ---------- + path_in is the path where the data tobe read are located + path_out is the path of teh directory where resulting files must be written + num output number + name a string which is used to write the names of the various files + thres_dens density threshold above which cells are considered + thres_level level threshold above which cells are considered + pos_zoom the center of the zoom coordinates + size_zoom the 3 zoom extension (x, y and z) + """ + + if path_out is not None: + directory_out = path_out + else: + directory_out = path_in + + name_txt = name + ".txt" + + # check whether hop entry files have been created (not test is done on .txt only + if len(glob.glob(name_txt)) == 0 or force: + ro = pymses.RamsesOutput(path_in, num) + amr = ro.amr_source(["rho"], grav_compat=gcomp) # density only is used + center = np.asarray(pos_zoom) + radius = size_zoom + min_coords = np.zeros(3) + max_coords = np.zeros(3) + min_coords[:] = center[:] - radius / 2.0 + max_coords[:] = center[:] + radius / 2.0 + + region = Box((min_coords, max_coords)) + + # region = Sphere(center,radius) + filt_amr = RegionFilter(region, amr) + cell_source = CellsToPoints( + filt_amr, + ) + + # selection of the cells of interest + def cell_selec_func(dset): + mask1 = dset["rho"] >= thres_dens + dx = dset.get_sizes() + mask2 = dx <= 0.5 ** thres_level + return mask1 * mask2 + + # begin cell_selec + cells_selec = PointFunctionFilter(cell_selec_func, cell_source).flatten() + dx = cells_selec.get_sizes() + ncells = cells_selec.npoints + + # fill the matrice with ID, x,y,z and masses of particles + val_mat = np.zeros((ncells, 5)) + + val_mat[:, 0] = np.arange(ncells) + val_mat[:, 1:4] = cells_selec.points[:, 0:3] + val_mat[:, 4] = cells_selec["rho"] * (dx ** 3) + + # write name.txt + head = str(ncells) + + np.savetxt( + name_txt, + val_mat, + fmt="%10d %.10e %.10e %.10e %.10e", + header=head, + delimiter=" ", + comments=" ", + ) + # end of creation name.txt + + # creation name.den + f = open(name + ".den", "wb") + f.write(pack("I", ncells)) + cells_selec["rho"].astype("f").tofile(f) + f.close() + + print(name + ".den created") + # end of creation name.den + + # HOP Algorithm + print("creation of .hop and .gbound du to hop") + fname = path_hop + name + ".txt" + print("look for hop in ", fname) + + h = HOP(fname, path_hop) + h.process_hop() + print("hop grouping is finished") + # end of HOP Algorithm + + idpart = val_mat[:, 0] + X = val_mat[:, 1] + Y = val_mat[:, 2] + Z = val_mat[:, 3] + mass = val_mat[:, 4] + + # read the gbound file to get list of particle numbers within groups + f = open(name + ".gbound", "r") + aline = f.readline() + ngroups = int(aline) + npart_v = np.zeros(ngroups, dtype=int) + for i in range(10): + aline = f.readline() + for i in range(ngroups): + aline = f.readline() + vec = aline.split() + igroup = int(vec[0]) + npart_v[igroup] = int(vec[1]) + f.close() + + # get the igroup array + group_ids = h.get_group_ids() + + # sort it and apply the sorting to the coordinates + # this means that the particules of group 1 are written first then of group 2 etc... + ind_sort = np.argsort(group_ids) + xx_v = X[ind_sort] + yy_v = Y[ind_sort] + zz_v = Z[ind_sort] + vect_id_group = group_ids[ind_sort] # not so useful + + # write the sorted cells + name_save_clump = directory_out + name + ".save" + + np.savez( + name_save_clump, + ngroups=ngroups, + npart_v=npart_v, + xx_v=xx_v, + yy_v=yy_v, + zz_v=zz_v, + vect_id_group=vect_id_group, + num=num, + name=name, + thres_dens=thres_dens, + thres_level=thres_level, + pos_zoom=pos_zoom, + size_zoom=size_zoom, + ) + + return name_save_clump diff --git a/pipeline_disk.py b/pipeline_disk.py index fe402c3..68b20b6 100644 --- a/pipeline_disk.py +++ b/pipeline_disk.py @@ -57,17 +57,36 @@ parser.add_argument( 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 ) +parser.add_argument( + "--binning", + help="Kind of binning (logarithmic or linear)", + choices=["log", "lin"], + default="log", +) +parser.add_argument( + "--rad_ext", help="Value of the highest bin", type=float, default=1.0 +) parser.add_argument( "-x", help="x position of the central point", type=float, default=1.0 ) @@ -125,13 +144,15 @@ for run in runs: while not success: try: + maps_disk = None if args.maps: - dp.make_image_disk( + print("[{}, {}] computing maps".format(run, i)) + maps_disk = dp.make_image_disk( path_in, i, path_out=path_out, tag=run, - map_size=1024, + map_size=args.mapsize, force=args.force_redo, level=args.level, cpu=args.cpu, @@ -139,17 +160,22 @@ for run in runs: fft=args.fft, interactive=args.interactive, ) - if args.disk or args.compare or args.evolution: + print("[{}, {}] maps computed".format(run, i)) + if args.disk: + print("[{}, {}] computing disk props".format(run, i)) dp.disk_prop( path_in, i, path_out=path_out, - rad_ext=1, nb_bin=args.nb_bin, + binning=args.binning, + rad_ext=args.rad_ext, force=args.force_redo, pos_star=np.array([args.x, args.y, args.z]), ) + print("[{}, {}] disk_props computed".format(run, i)) if args.disk: + print("[{}, {}] plotting disk props".format(run, i)) dp.plot_disk_prop( path_out, i, @@ -157,14 +183,26 @@ for run in runs: force=args.force_redo, interactive=args.interactive, ) + print("[{}, {}] disk props plotted".format(run, i)) + if args.pdf: + print("[{}, {}] computing pdf".format(run, i)) + dp.disk_pdf( + path_out, + i, + maps_disk, + pos_star=np.array([args.x, args.y, args.z]), + force=args.force_redo, + tag=run, + interactive=args.interactive, + ) + print("[{}, {}] pdf computed".format(run, i)) success = True - except ValueError: + except ValueError as e: + print(e) if args.watch and failures < args.allowed_failures: failures = failures + 1 print( - "Unable to proceed for run {} \ - output {}. Trying again in {} s ({} \ - tries remaining)".format( + "Unable to proceed for run {} output {}. Trying again in {} s ({} tries remaining)".format( run, i, args.waiting_time, args.allowed_failures - failures ) )