diff --git a/disk_postprocess.py b/disk_postprocess.py index 9387626..e151076 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -887,6 +887,15 @@ def disk_pdf( 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" @@ -1007,19 +1016,32 @@ def disk_pdf( P.close() -def compare(path, runs, num, force=False, interactive=False, Q_in_name=True, pdf=False): +def compare( + path, + runs, + nums, + force=False, + interactive=False, + Q_in_name=True, + pdf=False, + gamma=5.0 / 3.0, +): """ - Compare properties of a disk in several simulations + Compare time averaged properties of a disk in several simulations - num id of the ramses output + 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 - time = np.zeros(len(runs)) alpha_rey = np.zeros(len(runs)) alpha_grav = np.zeros(len(runs)) Q = np.zeros(len(runs)) @@ -1031,34 +1053,34 @@ def compare(path, runs, num, force=False, interactive=False, Q_in_name=True, pdf for i, run in enumerate(runs): path_run = path + "/" + run - # 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() + 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) - time[i] = prop_disk["time"] - alpha_rey[i] = prop_disk["alpha_rey_mean"] - alpha_grav[i] = prop_disk["alpha_grav_mean"] - Q[i] = prop_disk["Q_min"] if Q_in_name: Q0[i] = float(run.split("_")[2][1:]) - if pdf: - fit = prop_disk["fit"] - beta[i] = fit["beta"] - slope[i] = fit["slope"] # Check if the output file exists, and exit if it is the case - name_save = path + "/alphaQ_" + str(num).zfill(5) + out_ext + name_save = path + "/alphaQ_" + nums_name + out_ext # if (not force and len(glob.glob(name_save)) !=0): # return - title = "t=" + str(time[0]) + " (code)" - # alpha = f(Qmin) P.yscale("log") P.ylim([1e-7, 1.0]) @@ -1074,7 +1096,7 @@ def compare(path, runs, num, force=False, interactive=False, Q_in_name=True, pdf if interactive: P.figure() else: - P.savefig(path + "/alphaQ_" + str(num).zfill(5) + out_ext) + P.savefig(path + "/alphaQ_" + nums_name + out_ext) P.close() if Q_in_name: @@ -1093,7 +1115,7 @@ def compare(path, runs, num, force=False, interactive=False, Q_in_name=True, pdf if interactive: P.figure() else: - P.savefig(path + "/alphaQ0_" + str(num).zfill(5) + out_ext) + P.savefig(path + "/alphaQ0_" + nums_name + out_ext) P.close() if pdf: @@ -1108,7 +1130,7 @@ def compare(path, runs, num, force=False, interactive=False, Q_in_name=True, pdf if interactive: P.figure() else: - P.savefig(path + "/dcolslopebeta_" + str(num).zfill(5) + out_ext) + P.savefig(path + "/dcolslopebeta_" + nums_name + out_ext) P.close() # alpha = f(beta) @@ -1116,8 +1138,12 @@ def compare(path, runs, num, force=False, interactive=False, Q_in_name=True, pdf 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}$") @@ -1126,11 +1152,11 @@ def compare(path, runs, num, force=False, interactive=False, Q_in_name=True, pdf if interactive: P.figure() else: - P.savefig(path + "/alphabeta_" + str(num).zfill(5) + out_ext) + P.savefig(path + "/alphabeta_" + nums_name + out_ext) P.close() -def evolution(path, nums, force=False, interactive=False): +def evolution(path, nums, force=False, interactive=False, pdf=False): """ Plot properties over time @@ -1148,6 +1174,8 @@ def evolution(path, nums, force=False, interactive=False): 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): @@ -1168,6 +1196,8 @@ def evolution(path, nums, force=False, interactive=False): 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 @@ -1222,161 +1252,15 @@ def evolution(path, nums, force=False, interactive=False): 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 + # 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() diff --git a/pipeline_disk.py b/pipeline_disk.py index 3f241b2..3ee93f2 100644 --- a/pipeline_disk.py +++ b/pipeline_disk.py @@ -6,7 +6,7 @@ import argparse import time import numpy as np import disk_postprocess as dp - +from functools import reduce storage_in = "/home/nbrucy/simus/" storage_out = "/home/nbrucy/visus/" @@ -97,6 +97,10 @@ parser.add_argument( "-z", help="z position of the central point", type=float, default=1.0 ) +parser.add_argument( + "-ta", "--time_avg", help="Plot time averaged comparaison", action="store_true" +) + parser.add_argument( "--colormap", help="Colormap used", choices=dp.P.colormaps(), default="plasma" @@ -125,6 +129,9 @@ dp.out_ext = "." + args.format dp.P.rcParams["image.cmap"] = args.colormap dp.P.rcParams["savefig.dpi"] = args.dpi +# List of id that were successfully computed +run_succeded = {} + for run in runs: path_suffix = project + "/" + run path_in = storage_in + path_suffix @@ -138,6 +145,8 @@ for run in runs: except: pass + run_succeded[run] = [] + for i in range(first, last + 1, step): failures = 0 success = False @@ -196,8 +205,10 @@ for run in runs: interactive=args.interactive, ) print("[{}, {}] pdf computed".format(run, i)) + # If we are here, success ! success = True - except ValueError as e: + run_succeded[run].append(i) + except (ValueError, IOError) as e: print(e) if args.watch and failures < args.allowed_failures: failures = failures + 1 @@ -212,24 +223,45 @@ for run in runs: else: raise if args.evolution: + print("[{}] plotting evolution".format(run)) dp.evolution( path_out, - range(first, last + 1, step), + run_succeded[run], force=args.force_redo, interactive=args.interactive, + pdf=args.pdf, ) + print("[{}] evolution plotted".format(run)) if args.compare: path_suffix = project path = storage_out + path_suffix - for i in range(first, last + 1, step): - dp.compare( - path, - runs, - i, - force=args.force_redo, - interactive=args.interactive, - Q_in_name=(not args.pdf), - pdf=args.pdf, - ) + if args.time_avg: + # Select output availables for all runs + output_ok = reduce(np.intersect1d, [run_succeded[run] for run in runs]) + print(output_ok) + if len(output_ok) >= 1: + dp.compare( + path, + runs, + output_ok, + force=args.force_redo, + interactive=args.interactive, + Q_in_name=(not args.pdf), + pdf=args.pdf, + ) + else: + for i in range(first, last + 1, step): + # Select usable runs + runs_ok = [run for run in runs if i in run_succeded[run]] + if len(runs_ok) > 1: + dp.compare( + path, + runs_ok, + i, + force=args.force_redo, + interactive=args.interactive, + Q_in_name=(not args.pdf), + pdf=args.pdf, + )