diff --git a/disk_postprocess.py b/disk_postprocess.py index ef99e73..3c97d95 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -3,121 +3,161 @@ import sys import numpy as np import os import pymses + +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 import pickle as pickle -import module_extract as me -from pymses.filters import CellsToPoints from pymses.sources.ramses import output from pymses.analysis import Camera, raytracing, slicing +from pymses.filters import CellsToPoints +from pymses.analysis import ScalarOperator, FractionOperator, MaxLevelOperator + +P.rcParams["image.cmap"] = "plasma" +P.rcParams["savefig.dpi"] = 800 -def make_image_temp( - amr, - ro, - center, - radius, - num, +def make_image_disk( path, - i_im=0, - force=False, + num, path_out=None, - col_dens_only=False, - map_size=512, - vel_red=20, - im_pos=0.0, + order="<", + force=False, tag="", - cpuamr=False, + vel_red=20, + map_size=512, put_title=True, - mass_norm=1.0, - AU_units=False, + cpuamr=False, ): """ Make several useful image of an output of a simulation Parameters ---------- - amr - ro pymses.RamsesOutput object - - center 3D array for coordinates center - num output number + path path of the Ramses output + num Ramses output number + path_out path of the pipeline output + order '<' or '>' TODO + 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 """ - lbox = ro.info["boxlen"] # boxlen in codeunits (=>pc) + ro = pymses.RamsesOutput(path, num, order=order) + amr = ro.amr_source(["rho", "vel", "P", "Bl", "Br"]) + rad = 0.5 + center = [0.5, 0.5, 0.5] - ( - AU, - pc, - Ms, - Myr, - scale_n, - scale_d, - scale_t, - scale_l, - scale_v, - scale_T2, - scale_ener, - scale_mag, - microG, - km_s, - Cwnm, - scale_mass, - unit_col, - lbox_pc, - ) = normalisation(ro) + make_image_aux( + amr, + ro, + center, + rad, + num, + path, + force=force, + path_out=path_out, + map_size=map_size, + vel_red=vel_red, + tag=tag, + cpuamr=cpuamr, + put_title=put_title, + ) + + +def make_image_aux( + amr, + ro, + center, + radius, + num, + path, + force=False, + path_out=None, + map_size=512, + vel_red=20, + tag="", + cpuamr=False, + put_title=True, +): + """ + 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 + cpuamr plot also levels and cpus at each step + """ + + lbox = ro.info["boxlen"] # boxlen in codeunits + lbox_units = lbox + + G = 1.0 # Gravitational constant ntick = 6 - if AU_units: - units = pc / AU - lbox_units = lbox * units - titre_x = "x (AU)" - titre_y = "y (AU)" - titre_z = "z (AU)" - else: - units = 1.0 - lbox_units = lbox - titre_x = "x (pc)" - titre_y = "y (pc)" - titre_z = "z (pc)" - - lbox_cm = lbox * pc # lbox in cm + title_ax = {"x": "x (code)", "y": "y (code)", "z": "z (code)"} time = ro.info["time"] # time in codeunits - time = time * scale_t / Myr - - titre = "t=" + str(time)[0:5] + " (Myr)" + title = "t=" + str(time)[0:5] + " (code)" if path_out is not None: directory = path_out else: directory = path - rt = raytracing.RayTracer(amr, ro.info, rho_op) - rho_op = ScalarOperator(lambda dset: dset["rho"], ro.info["unit_density"]) + name = directory + "/coldens_z" + "_" + tag + "_" + format(num, "05") + ".jpeg" + if len(glob.glob(name)) == 1 and not force: + return - ax_names = ["x", "y", "z"] - up_axes = ["z", "z", "x"] - other_axes = ["z", "x", "y"] + rho_op = ScalarOperator(lambda dset: dset["rho"], ro.info["unit_density"]) + 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} image_names = ["coldens", "rho", "T"] - for i, ax_name in enumerate(ax_names): + 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_name, + line_of_sight_axis=ax_los, region_size=[2.0 * radius, 2.0 * radius], distance=radius, far_cut_depth=radius, - up_vector=up_axes[i], + up_vector=ax_v, map_max_size=map_size, ) datamap = rt.process(cam, surf_qty=True) - map_col = np.log10(datamap.map.T * lbox_cm) - P.clf() + # Column density + dmap_col = datamap.map.T * lbox + map_col = np.log10(dmap_col) + + P.close() im = P.imshow( map_col, extent=[ @@ -129,38 +169,44 @@ def make_image_temp( origin="lower", ) - P.locator_params(axis=other_axes[i], nbins=ntick) - P.locator_params(axis=up_axes[i], nbins=ntick) + P.locator_params(axis=ax_h, nbins=ntick) + P.locator_params(axis=ax_v, nbins=ntick) if put_title: - P.title(titre) + P.title(title) + + P.xlabel(title_ax[ax_h]) + P.ylabel(title_ax[ax_v]) - P.xlabel(titre_x) - P.ylabel(titre_y) cbar = P.colorbar(im) - cbar.set_label(r"$log(N) \, cm^{-2}$") - - name = directory + "/coldens_" + ax_name + "_" + tag + +"_" + format(num, "05") + cbar.set_label(r"$log(N)$ (code)") + name = directory + "/coldens_" + ax_los + "_" + tag + "_" + format(num, "05") name_im = name + ".jpeg" P.savefig(name_im) + # 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 - vx_op = ScalarOperator(lambda dset: dset["vel"][:, 0], ro.info["unit_velocity"]) - dmap_vx = slicing.SliceMap(amr, cam_z, vx_op, z=0.0) - map_vx_red = dmap_vx.map[::vel_red, ::vel_red] + 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_vx_red = map_vx_red.T + map_vh_red = map_vh_red.T - vy_op = ScalarOperator(lambda dset: dset["vel"][:, 1], ro.info["unit_velocity"]) - dmap_vy = slicing.SliceMap(amr, cam, vy_op, z=0.0) - map_vy_red = dmap_vy.map[::vel_red, ::vel_red] + 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 - map_vy_red = map_vy_red.T - - P.clf() + P.close() im = P.imshow( map_rho, extent=[ @@ -171,50 +217,45 @@ def make_image_temp( ], origin="lower", ) - - P.locator_params(axis=other_axes[i], nbins=ntick) - P.locator_params(axis=up_axes[i], nbins=ntick) - - nx = map_vx_red.shape[0] - ny = map_vx_red.shape[1] - vec_x = ( - np.arange(nx) * 2.0 / nx * radius - radius + center[0] + radius / nx + 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_units - vec_y = ( - np.arange(ny) * 2.0 / ny * radius - radius + center[1] + radius / nx + vec_v = ( + np.arange(nv) * 2.0 / nv * radius - radius + center[1] + radius / nv ) * lbox_units - xx, yy = np.meshgrid(vec_x, vec_y) - - map_vx_red = map_vx_red * scale_v / km_s - map_vy_red = map_vy_red * scale_v / km_s - max_v = np.max(np.sqrt(map_vx_red ** 2 + map_vy_red ** 2)) - - Q = P.quiver(xx, yy, map_vx_red, map_vy_red, units="width") + 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] + "\, km \, s^{-1}$", + r"$" + str(max_v)[0:4] + "$ (code)", labelpos="E", coordinates="figure", ) if put_title: - P.title(titre) - P.xlabel(titre_x) - P.ylabel(titre_y) + 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) \, (cm^{-3})$") + cbar.set_label(r"$log(n)$ (code)") - name = directory + "/rho_z" + "_" + tag + str(i_im) + "_" + format(num, "05") + name = directory + "/rho_" + ax_los + "_" + tag + "_" + format(num, "05") name_im = name + ".jpeg" P.savefig(name_im) P_op = ScalarOperator(lambda dset: dset["P"], ro.info["unit_pressure"]) - dmap_P = slicing.SliceMap(amr, cam_z, P_op, z=0.0) - P.clf() - map_T = np.log10(dmap_P.map.T / dmap_rho.map.T * scale_T2) + 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) im = P.imshow( map_T, @@ -231,25 +272,56 @@ def make_image_temp( P.locator_params(axis="y", nbins=ntick) if put_title: - P.title(titre) - - P.xlabel(titre_x) - P.ylabel(titre_y) + 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_z" + "_" + tag + str(i_im) + "_" + format(num, "05") + name = directory + "/T_" + ax_los + "_" + tag + "_" + format(num, "05") name_im = name + ".jpeg" 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 + ) + + im = P.imshow( + map_T, + extent=[ + (-radius + center[0]) * lbox_units, + (radius + center[0]) * lbox_units, + (-radius + center[1]) * lbox_units, + (radius + center[1]) * lbox_units, + ], + 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 + "/Q_" + ax_los + "_" + tag + "_" + format(num, "05") + name_im = name + ".jpeg" + P.savefig(name_im) + if cpuamr: level_op = MaxLevelOperator() amr.set_read_levelmax(20) rt = raytracing.RayTracer(amr, ro.info, level_op) - datamap = rt.process(cam_z, surf_qty=True) + datamap = rt.process(cam, surf_qty=True) map_level = datamap.map.T - P.clf() + P.close() im = P.imshow( map_level, extent=[ @@ -265,49 +337,25 @@ def make_image_temp( P.locator_params(axis="y", nbins=ntick) if put_title: - P.title(titre) - P.xlabel(titre_x) - P.ylabel(titre_y) + P.title(title) + P.xlabel(title_ax[ax_h]) + P.ylabel(title_ax[ax_v]) cbar = P.colorbar(im) cbar.set_label(r"level") - name = ( - directory + "/level_z" + "_" + tag + str(i_im) + "_" + format(num, "05") - ) + name = directory + "/level_" + ax_los + "_" + tag + "_" + format(num, "05") name_im = name + ".jpeg" P.savefig(name_im) - if ps: - P.savefig(name + ".ps") - - if descrip is not None: - dd = {} - dd.update( - { - "name": "mean AMR level in the z-direction" - + " (" - + tag - + str(i_im) - + ")" - } - ) - dd.update({"type": "image"}) - dd.update({"description": "Mean AMR level in the z-direction."}) - dd.update({"display": 0}) - dd.update({"item": i_im}) - dd.update({"tag": tag}) - data = {"jpeg": name_im} - dd.update({"data": data}) - list_im.append(dd) 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_z, surf_qty=True) + datamap = rt.process(cam, surf_qty=True) map_cpu = datamap.map.T - P.clf() + P.close() im = P.imshow( map_cpu, extent=[ @@ -322,16 +370,14 @@ def make_image_temp( P.locator_params(axis="y", nbins=ntick) if put_title: - P.title(titre) + P.title(title) - P.xlabel(titre_x) - P.ylabel(titre_y) + P.xlabel(title_ax[ax_h]) + P.ylabel(title_ax[ax_v]) cbar = P.colorbar(im) - cbar.set_label(r"level") + cbar.set_label(r"cpu") - name = ( - directory + "/cpu_z" + "_" + tag + str(i_im) + "_" + format(num, "05") - ) + name = directory + "/cpu_" + ax_los + "_" + tag + "_" + format(num, "05") name_im = name + ".jpeg" P.savefig(name_im) @@ -409,8 +455,7 @@ def disk_prop( v_az = (pos[:, 0] * vel[:, 1] - pos[:, 1] * vel[:, 0]) / norm_pos # Select cells that are actually in the disk, ie within the scale height - # TODO Check units - G = 1.0 # G=6.8e-8 + 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 @@ -439,6 +484,7 @@ def disk_prop( v_kepl_rad = np.zeros(nb_bin - 1) v_rad_rad = np.zeros(nb_bin - 1) alpha_rey_rad = np.zeros(nb_bin - 1) + Q_kepl_rad = np.zeros(nb_bin - 1) for i in range(nb_bin - 1): mask_bin = (rc_disk > rad[i]) & (rc_disk < rad[i + 1]) @@ -502,9 +548,11 @@ def disk_prop( v_rad_rad = v_rad_rad # * scale_v / km_s v_kepl_rad = v_kepl_rad + Q_kepl_rad = cs_rad * v_az_rad / (np.pi * G * coldens_rad * rad[0 : nb_bin - 1]) + prop_disk = { "time": time, - "rad_AU": rad[0 : nb_bin - 1], + "rad": rad[0 : nb_bin - 1], "center": pos_star, "alpha_rey": alpha_rey_rad, "v_rad": v_rad_rad, @@ -515,6 +563,7 @@ def disk_prop( "press": press_rad, "temp": temp_rad, "cs": cs_rad, + "Q_kepl": Q_kepl_rad, } # store the results @@ -545,12 +594,14 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): prop_disk = pickle.load(f) f.close() - P.clf() + # Check if the output file exists, and exit if it is the case + name_save = path + "/rho_disk_r_" + str(num).zfill(5) + ".jpeg" + if not force and len(glob.glob(name_save)) != 0: + return + + P.close() P.plot( - np.log10(prop_disk["rad_AU"]), - np.log10(prop_disk["rho"]), - color="k", - linewidth=2, + np.log10(prop_disk["rad"]), np.log10(prop_disk["rho"]), color="k", linewidth=2 ) P.ylabel(r"$\log(n) \, (cm^{-3})$") P.xlabel("disk radius") @@ -559,12 +610,9 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + ".pdf") P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + ".jpeg") - P.clf() + P.close() P.plot( - np.log10(prop_disk["rad_AU"]), - np.log10(prop_disk["temp"]), - color="k", - linewidth=2, + np.log10(prop_disk["rad"]), np.log10(prop_disk["temp"]), color="k", linewidth=2 ) P.ylabel(r"$\log(T) \, (K)$") P.xlabel("disk radius") @@ -573,15 +621,15 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + ".pdf") P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + ".jpeg") - P.clf() + P.close() P.xscale("log") P.yscale("symlog", linthreshy=0.01) - P.plot((prop_disk["rad_AU"]), ((prop_disk["v_rad"])), color="k", linewidth=2) - P.plot((prop_disk["rad_AU"]), ((prop_disk["v_kepl"])), color="b", linewidth=2) - P.plot((prop_disk["rad_AU"]), (abs(prop_disk["v_az"])), color="r", linewidth=2) - P.plot((prop_disk["rad_AU"]), ((prop_disk["cs"])), color="c", linewidth=2) + P.plot((prop_disk["rad"]), ((prop_disk["v_rad"])), color="k", linewidth=2) + P.plot((prop_disk["rad"]), ((prop_disk["v_kepl"])), color="b", linewidth=2) + P.plot((prop_disk["rad"]), (abs(prop_disk["v_az"])), color="r", linewidth=2) + P.plot((prop_disk["rad"]), ((prop_disk["cs"])), color="c", linewidth=2) P.legend((r"$v_r$", r"$v_{kepl}$", r"$v_\phi$", r"$c_s$"), loc="upper right") @@ -592,9 +640,9 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + ".pdf") P.savefig(path + "/V_disk_r_" + str(num).zfill(5) + ".jpeg") - P.clf() + P.close() P.plot( - np.log10(prop_disk["rad_AU"]), + np.log10(prop_disk["rad"]), np.log10(prop_disk["coldens"]), color="k", linewidth=2, @@ -606,13 +654,13 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + ".pdf") P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + ".jpeg") - P.clf() + # Alpha + P.close() P.xscale("log") - P.yscale("symlog", linthreshy=0.001) + P.yscale("log") + P.ylim([1e-5, 1.0]) - P.plot(prop_disk["rad_AU"], abs(prop_disk["alpha_rey"]), color="b", linewidth=2) - - # P.legend(r'$\alpha_{Rey}$', loc='upper right') + P.plot(prop_disk["rad"], abs(prop_disk["alpha_rey"]), color="b", linewidth=2) P.ylabel(r"$\alpha$") P.xlabel("disk radius ") @@ -621,6 +669,17 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + ".pdf") P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + ".jpeg") + # Q + + P.close() + 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 ") + 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.savefig(path + "/Q_r_" + str(num).zfill(5) + ".pdf") + P.savefig(path + "/Q_r_" + str(num).zfill(5) + ".jpeg") diff --git a/pipeline_disk.py b/pipeline_disk.py index 058087b..62513c0 100644 --- a/pipeline_disk.py +++ b/pipeline_disk.py @@ -1,11 +1,11 @@ # coding: utf-8 import sys import numpy as np -import module_extract as me -import extract_disk as d import pymses import os +from shutil import copy import argparse +import time import disk_postprocess as dp @@ -24,19 +24,39 @@ 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("-p", "--project", help="specify project name", default="disk") + +parser.add_argument( + "-w", "--watch", help="wait and watch for missing outputs", action="store_true" +) +parser.add_argument( + "-wt", + "--waiting_time", + help="time between to successive try when watching new outputs (in second)", + type=int, + default=120, +) +parser.add_argument( + "-af", + "--allowed_failures", + help="number of allowed failures when waiting", + default=30, +) + args = parser.parse_args() user = "nbrucy" folder = "simus" -project = "disk" +project = args.project runs = args.runs first = args.first_output last = args.last_output step = args.step + for run in runs: path_suffix = user + "/" + folder + "/" + project + "/" + run path_in = storage_in + path_suffix @@ -44,23 +64,32 @@ for run in runs: if not os.path.exists(path_out): os.makedirs(path_out) + copy(path_in + "/disk.nml", path_out) + copy(path_in + "/output_00001/compilation.txt", path_out) for i in range(first, last + 1, step): - me.make_image_zoom( - path_in, - i, - [0.5], - sinks=False, - force=False, - path_out=path_out, - tag=run + "_", - cpuamr=False, - mag_im=False, - AU_units=False, - ) - # me.look(path_in, i) - if args.disk: - dp.disk_prop( - path_in, i, path_out=path_out, rad_ext=1, nb_bin=50, force=True - ) - dp.plot_disk_prop(path_out, i, tag=run + "_") + failures = 0 + success = False + + while not success: + try: + dp.make_image_disk( + path_in, i, path_out=path_out, tag=run, map_size=1024 + ) + if args.disk: + dp.disk_prop( + path_in, i, path_out=path_out, rad_ext=1, nb_bin=50, force=False + ) + dp.plot_disk_prop(path_out, i, tag=run) + success = True + except ValueError: + 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( + run, i, args.waiting_time, args.allowed_failures - failures + ) + ) + time.sleep(args.waiting_time) + else: + raise