diff --git a/disk_postprocess.py b/disk_postprocess.py index 3c97d95..31157ef 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -19,8 +19,11 @@ from pymses.analysis import Camera, raytracing, slicing 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"] = 800 +P.rcParams["savefig.dpi"] = 400 def make_image_disk( @@ -52,7 +55,7 @@ def make_image_disk( """ ro = pymses.RamsesOutput(path, num, order=order) - amr = ro.amr_source(["rho", "vel", "P", "Bl", "Br"]) + amr = ro.amr_source(["rho", "vel", "P"]) rad = 0.5 center = [0.5, 0.5, 0.5] @@ -86,6 +89,7 @@ def make_image_aux( vel_red=20, tag="", cpuamr=False, + pos_star=np.array([1.0, 1.0, 1.0]), put_title=True, ): """ @@ -122,7 +126,7 @@ def make_image_aux( else: directory = path - name = directory + "/coldens_z" + "_" + tag + "_" + format(num, "05") + ".jpeg" + name = directory + "/coldens_z" + "_" + tag + "_" + format(num, "05") + out_ext if len(glob.glob(name)) == 1 and not force: return @@ -181,7 +185,7 @@ def make_image_aux( cbar = P.colorbar(im) cbar.set_label(r"$log(N)$ (code)") name = directory + "/coldens_" + ax_los + "_" + tag + "_" + format(num, "05") - name_im = name + ".jpeg" + name_im = name + out_ext P.savefig(name_im) # Rho slice @@ -248,7 +252,7 @@ def make_image_aux( cbar.set_label(r"$log(n)$ (code)") name = directory + "/rho_" + ax_los + "_" + tag + "_" + format(num, "05") - name_im = name + ".jpeg" + name_im = name + out_ext P.savefig(name_im) P_op = ScalarOperator(lambda dset: dset["P"], ro.info["unit_pressure"]) @@ -279,18 +283,48 @@ def make_image_aux( cbar.set_label(r"$log(T) \, (K)$") name = directory + "/T_" + ax_los + "_" + tag + "_" + format(num, "05") - name_im = name + ".jpeg" + name_im = name + out_ext P.savefig(name_im) # Toomre parameter - if False and ax_los == "z": - map_Q = np.log10(np.sqrt(dmap_P.map.T / dmap_rho.map.T)) / ( - np.pi * map_col * G - ) + 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) + rt_cs = raytracing.RayTracer(amr, ro.info, cs_op) + + 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 = np.log10(dmap_Q) + + P.close() im = P.imshow( - map_T, + map_Q, extent=[ (-radius + center[0]) * lbox_units, (radius + center[0]) * lbox_units, @@ -305,14 +339,14 @@ def make_image_aux( 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)$") + P.xlabel(title_ax[ax_h]) + P.ylabel(title_ax[ax_v]) + cbar = P.colorbar(im) + cbar.set_label(r"$log(Q)$") - name = directory + "/Q_" + ax_los + "_" + tag + "_" + format(num, "05") - name_im = name + ".jpeg" - P.savefig(name_im) + name = directory + "/Q_" + ax_los + "_" + tag + "_" + format(num, "05") + name_im = name + out_ext + P.savefig(name_im) if cpuamr: level_op = MaxLevelOperator() @@ -344,7 +378,7 @@ def make_image_aux( cbar.set_label(r"level") name = directory + "/level_" + ax_los + "_" + tag + "_" + format(num, "05") - name_im = name + ".jpeg" + name_im = name + out_ext P.savefig(name_im) cpu_op = ScalarOperator( @@ -378,7 +412,7 @@ def make_image_aux( cbar.set_label(r"cpu") name = directory + "/cpu_" + ax_los + "_" + tag + "_" + format(num, "05") - name_im = name + ".jpeg" + name_im = name + out_ext P.savefig(name_im) @@ -432,7 +466,7 @@ def disk_prop( ro = pymses.RamsesOutput(path_in, num) lbox = ro.info["boxlen"] # boxlen in codeunits (=>pc) - time = ro.info["time"] # * scale_t / Myr + time = ro.info["time"] # time in codeunits # Get array of cell positions amr = ro.amr_source(["rho", "vel", "Br", "Bl", "P"]) @@ -440,9 +474,9 @@ def disk_prop( cells = cell_source.flatten() dx = cells.get_sizes() 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 @@ -474,6 +508,8 @@ def disk_prop( v_az_disk = v_az[mask] v_kepl = np.sqrt(mass_star * G / rc_disk) + total_mass_disk = np.sum(rho_disk * dvol_disk) + # Initialize binned quantities cs_rad = np.zeros(nb_bin - 1) temp_rad = np.zeros(nb_bin - 1) @@ -552,6 +588,7 @@ def disk_prop( prop_disk = { "time": time, + "tot_mass": total_mass_disk, "rad": rad[0 : nb_bin - 1], "center": pos_star, "alpha_rey": alpha_rey_rad, @@ -572,16 +609,13 @@ def disk_prop( f.close() -def plot_disk_prop(path, num, force=False, pdf=False, tag=""): +def plot_disk_prop(path, num, force=False, tag=""): """ 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 - - pdf if set, do output in pdf as well - """ # Load property file @@ -595,16 +629,22 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): 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) + ".jpeg" + 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["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 ) - P.ylabel(r"$\log(n) \, (cm^{-3})$") + 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 pdf: P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + ".pdf") @@ -616,6 +656,8 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): ) 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") @@ -636,9 +678,7 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): P.ylabel(r"$V \, (km s^{-1})$") P.xlabel("disk radius") - if pdf: - P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + ".pdf") - P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + ".jpeg") + P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + out_ext) P.close() P.plot( @@ -649,10 +689,8 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): ) P.ylabel(r"$\log(N) \, (cm^{-2})$") P.xlabel("disk radius ") - - if pdf: - P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + ".pdf") - P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + ".jpeg") + P.title(title) + P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + out_ext) # Alpha P.close() @@ -664,10 +702,8 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): P.ylabel(r"$\alpha$") P.xlabel("disk radius ") - - if pdf: - P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + ".pdf") - P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + ".jpeg") + P.title(title) + P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + out_ext) # Q @@ -675,10 +711,10 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): P.xscale("log") P.yscale("log", linthreshy=0.001) - 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 pdf: P.savefig(path + "/Q_r_" + str(num).zfill(5) + ".pdf") diff --git a/pipeline_disk.py b/pipeline_disk.py index 62513c0..eb42c87 100644 --- a/pipeline_disk.py +++ b/pipeline_disk.py @@ -14,18 +14,21 @@ storage_out = "/dsm/anais/storageA/" parser = argparse.ArgumentParser() parser.add_argument("runs", help="name of runs", nargs="*", default=["015_iso"]) -parser.add_argument( - "-f", "--first_output", help="id of first output", type=int, default=1 -) -parser.add_argument( - "-l", "--last_output", help="id of last output", type=int, default=100 -) +parser.add_argument("-b", "--begin", help="id of first output", type=int, default=1) +parser.add_argument("-e", "--end", help="id of last output", type=int, default=100) parser.add_argument("-s", "--step", help="step between two output", type=int, default=1) parser.add_argument( "-d", "--disk", help="do specific disk radial analysis", action="store_true" ) parser.add_argument("-p", "--project", help="specify project name", default="disk") +parser.add_argument( + "-fr", + "--force_redo", + help="redo plots even if the files already exist", + action="store_true", +) + parser.add_argument( "-w", "--watch", help="wait and watch for missing outputs", action="store_true" ) @@ -52,8 +55,8 @@ folder = "simus" project = args.project runs = args.runs -first = args.first_output -last = args.last_output +first = args.begin +last = args.end step = args.step @@ -74,13 +77,23 @@ for run in runs: while not success: try: dp.make_image_disk( - path_in, i, path_out=path_out, tag=run, map_size=1024 + path_in, + i, + path_out=path_out, + tag=run, + map_size=1024, + force=args.force_redo, ) if args.disk: dp.disk_prop( - path_in, i, path_out=path_out, rad_ext=1, nb_bin=50, force=False + path_in, + i, + path_out=path_out, + rad_ext=1, + nb_bin=50, + force=args.force_redo, ) - dp.plot_disk_prop(path_out, i, tag=run) + dp.plot_disk_prop(path_out, i, tag=run, force=args.force_redo) success = True except ValueError: if args.watch and failures < args.allowed_failures: