From faa8568bd5b8884971c0314b57fb887c7824aabc Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Thu, 18 Apr 2019 16:04:20 +0200 Subject: [PATCH] Correct vel/Q lbox factor, add interactive mode --- disk_postprocess.py | 161 ++++++++++++++++++++++++++++---------------- pipeline_disk.py | 32 ++++++--- 2 files changed, 126 insertions(+), 67 deletions(-) diff --git a/disk_postprocess.py b/disk_postprocess.py index 31157ef..da62519 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -37,6 +37,8 @@ def make_image_disk( map_size=512, put_title=True, cpuamr=False, + pos_star=np.array([1.0, 1.0, 1.0]), + interactive=False, ): """ Make several useful image of an output of a simulation @@ -45,7 +47,7 @@ def make_image_disk( ---------- path path of the Ramses output num Ramses output number - path_out path of the pipeline output + path_out path of the pipeline outputb order '<' or '>' TODO force if set, erase any existing pipeline output files tag string to add to the output name @@ -73,6 +75,8 @@ def make_image_disk( tag=tag, cpuamr=cpuamr, put_title=put_title, + pos_star=pos_star, + interactive=interactive, ) @@ -91,6 +95,7 @@ def make_image_aux( cpuamr=False, pos_star=np.array([1.0, 1.0, 1.0]), put_title=True, + interactive=False, ): """ Make several useful image of an output of a simulation, auxillary function @@ -161,7 +166,11 @@ def make_image_aux( dmap_col = datamap.map.T * lbox map_col = np.log10(dmap_col) - P.close() + if interactive: + P.figure() + else: + P.close() + im = P.imshow( map_col, extent=[ @@ -186,7 +195,12 @@ def make_image_aux( cbar.set_label(r"$log(N)$ (code)") name = directory + "/coldens_" + ax_los + "_" + tag + "_" + format(num, "05") name_im = name + out_ext - P.savefig(name_im) + + if interactive: + P.figure() + else: + P.savefig(name_im) + P.close() # Rho slice dmap_rho = slicing.SliceMap(amr, cam, rho_op, z=0.0) @@ -210,7 +224,6 @@ def make_image_aux( map_vv_red = dmap_vv.map[::vel_red, ::vel_red] map_vv_red = map_vv_red.T - P.close() im = P.imshow( map_rho, extent=[ @@ -253,11 +266,15 @@ def make_image_aux( name = directory + "/rho_" + ax_los + "_" + tag + "_" + format(num, "05") name_im = name + out_ext - P.savefig(name_im) + 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) - P.close() + dmap_T = dmap_P.map.T / dmap_rho.map.T map_T = np.log10(dmap_T) @@ -284,7 +301,11 @@ def make_image_aux( name = directory + "/T_" + ax_los + "_" + tag + "_" + format(num, "05") name_im = name + out_ext - P.savefig(name_im) + if interactive: + P.figure() + else: + P.savefig(name_im) + P.close() # Toomre parameter @@ -322,7 +343,6 @@ def make_image_aux( dmap_Q = (lbox * dmap_cs.map.T) * dmap_omega.map.T / (np.pi * G * dmap_col) map_Q = np.log10(dmap_Q) - P.close() im = P.imshow( map_Q, extent=[ @@ -346,16 +366,19 @@ def make_image_aux( name = directory + "/Q_" + ax_los + "_" + tag + "_" + format(num, "05") name_im = name + out_ext - P.savefig(name_im) + if interactive: + P.figure() + else: + P.savefig(name_im) + P.close() if cpuamr: level_op = MaxLevelOperator() amr.set_read_levelmax(20) - rt = raytracing.RayTracer(amr, ro.info, level_op) - datamap = rt.process(cam, surf_qty=True) + rt_level = raytracing.RayTracer(amr, ro.info, level_op) + datamap = rt_level.process(cam, surf_qty=True) map_level = datamap.map.T - P.close() im = P.imshow( map_level, extent=[ @@ -379,17 +402,21 @@ def make_image_aux( name = directory + "/level_" + ax_los + "_" + tag + "_" + format(num, "05") name_im = name + out_ext - P.savefig(name_im) + + if interactive: + P.figure() + else: + P.savefig(name_im) + P.close() cpu_op = ScalarOperator( lambda dset: dset.icpu * (np.ones(dset["P"].shape)), ro.info["unit_pressure"], ) - rt = raytracing.RayTracer(amr, ro.info, cpu_op) - datamap = rt.process(cam, surf_qty=True) + rt_cpu = raytracing.RayTracer(amr, ro.info, cpu_op) + datamap = rt_cpu.process(cam, surf_qty=True) map_cpu = datamap.map.T - P.close() im = P.imshow( map_cpu, extent=[ @@ -413,7 +440,12 @@ def make_image_aux( name = directory + "/cpu_" + ax_los + "_" + tag + "_" + format(num, "05") name_im = name + out_ext - P.savefig(name_im) + + if interactive: + P.figure() + else: + P.savefig(name_im) + P.close() def disk_prop( @@ -480,7 +512,7 @@ def disk_prop( # Get cylindrical radius rc = np.sqrt(pos[:, 0] ** 2 + pos[:, 1] ** 2) # Get velocities - vel = cells["vel"] * lbox + vel = cells["vel"] # Get radial component of velocity norm_pos = rc norm_pos[rc == 0] = 1.0e-10 # Avoid division per 0 @@ -537,7 +569,9 @@ def disk_prop( rho_rad[i] = np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) / np.sum( dvol_disk[mask_bin] ) - temp_rad[i] = press_rad[i] / rho_rad[i] + temp_rad[i] = np.sum(press_disk[mask_bin] * dvol_disk[mask_bin]) / np.sum( + rho_disk[mask_bin] * dvol_disk[mask_bin] + ) # TODO verifier unites # Surface of a bin : S = dr * 2 * pi * r with @@ -609,7 +643,7 @@ def disk_prop( f.close() -def plot_disk_prop(path, num, force=False, tag=""): +def plot_disk_prop(path, num, force=False, tag="", interactive=False): """ Plot properties of a disk @@ -623,7 +657,7 @@ def plot_disk_prop(path, num, force=False, tag=""): # Check if the properties file exists if len(glob.glob(name_save)) == 0: - throw("no pickle file for disk properties. Run single_z_disk_prop") + raise ("no pickle file for disk properties. Run single_z_disk_prop") f = open(name_save, "r") prop_disk = pickle.load(f) f.close() @@ -637,33 +671,32 @@ def plot_disk_prop(path, num, force=False, tag=""): mass = prop_disk["tot_mass"] title = "t=" + str(time)[0:5] + " (code)" - P.close() - P.plot( - np.log10(prop_disk["rad"]), np.log10(prop_disk["rho"]), color="k", linewidth=2 - ) + if interactive: + P.figure() + else: + P.close() + + P.xscale("log") + P.plot(prop_disk["rad"], np.log10(prop_disk["rho"]), color="k", linewidth=2) P.ylabel(r"$\log(n) \, (code)$") P.xlabel("disk radius") P.title(title) - P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + out_ext) + if interactive: + P.figure() + else: + P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + out_ext) + P.close() - if pdf: - P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + ".pdf") - P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + ".jpeg") - - P.close() - P.plot( - np.log10(prop_disk["rad"]), np.log10(prop_disk["temp"]), color="k", linewidth=2 - ) + P.xscale("log") + P.plot(prop_disk["rad"], np.log10(prop_disk["temp"]), color="k", linewidth=2) P.ylabel(r"$\log(T) \, (K)$") P.xlabel("disk radius") P.title(title) - P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + out_ext) - - if pdf: - P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + ".pdf") - P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + ".jpeg") - - P.close() + 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) @@ -677,44 +710,54 @@ def plot_disk_prop(path, num, force=False, tag=""): 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.savefig(path + "/V_disk_r_" + str(num).zfill(5) + out_ext) - - P.close() - P.plot( - np.log10(prop_disk["rad"]), - np.log10(prop_disk["coldens"]), - color="k", - linewidth=2, - ) + P.xscale("log") + P.plot(prop_disk["rad"], np.log10(prop_disk["coldens"]), color="k", linewidth=2) P.ylabel(r"$\log(N) \, (cm^{-2})$") P.xlabel("disk radius ") P.title(title) - P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + out_ext) + if interactive: + P.figure() + else: + P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + out_ext) + P.close() # Alpha - P.close() P.xscale("log") P.yscale("log") P.ylim([1e-5, 1.0]) P.plot(prop_disk["rad"], abs(prop_disk["alpha_rey"]), color="b", linewidth=2) + P.plot(prop_disk["rad"], abs(prop_disk["alpha_rey"]), color="b", linewidth=2) + P.ylabel(r"$\alpha$") P.xlabel("disk radius ") P.title(title) - P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + out_ext) + + if interactive: + P.figure() + else: + P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + out_ext) + P.close() # Q - - P.close() - P.xscale("log") - P.yscale("log", linthreshy=0.001) - + P.ylim([0, 10.0]) + P.plot(prop_disk["rad"], abs(prop_disk["Q_kepl"]), color="b", linewidth=2) P.ylabel(r"$Q$") P.xlabel("disk radius ") P.title(title + ", mass of disk = {} (code)".format(mass)) - P.savefig(path + "/Q_r_" + str(num).zfill(5) + out_ext) + + if interactive: + pass + else: + P.savefig(path + "/Q_r_" + str(num).zfill(5) + out_ext) + P.close() if pdf: P.savefig(path + "/Q_r_" + str(num).zfill(5) + ".pdf") diff --git a/pipeline_disk.py b/pipeline_disk.py index eb42c87..4d425ec 100644 --- a/pipeline_disk.py +++ b/pipeline_disk.py @@ -20,6 +20,8 @@ parser.add_argument("-s", "--step", help="step between two output", type=int, de 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("-p", "--project", help="specify project name", default="disk") parser.add_argument( @@ -45,6 +47,16 @@ parser.add_argument( help="number of allowed failures when waiting", default=30, ) +parser.add_argument("--cpuamr", help="plot levels and cpu", action="store_true") +parser.add_argument( + "-x", help="x position of the central point", type=float, default=1.0 +) +parser.add_argument( + "-y", help="y position of the central point", type=float, default=1.0 +) +parser.add_argument( + "-z", help="z position of the central point", type=float, default=1.0 +) args = parser.parse_args() @@ -76,14 +88,17 @@ for run in runs: while not success: try: - dp.make_image_disk( - path_in, - i, - path_out=path_out, - tag=run, - map_size=1024, - force=args.force_redo, - ) + if args.maps: + dp.make_image_disk( + path_in, + i, + path_out=path_out, + tag=run, + map_size=1024, + force=args.force_redo, + cpuamr=args.cpuamr, + pos_star=np.array([args.x, args.y, args.z]), + ) if args.disk: dp.disk_prop( path_in, @@ -92,6 +107,7 @@ for run in runs: rad_ext=1, nb_bin=50, force=args.force_redo, + pos_star=np.array([args.x, args.y, args.z]), ) dp.plot_disk_prop(path_out, i, tag=run, force=args.force_redo) success = True