diff --git a/disk_postprocess.py b/disk_postprocess.py index 3c961a5..0ba4cef 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -1,11 +1,13 @@ +# coding: utf-8 import sys import numpy as np import os import pymses import pylab as P import glob as glob -import pandas as pd import pickle as pickle +import module_extract as me +from pymses.filters import CellsToPoints from pymses.sources.ramses import output @@ -38,7 +40,7 @@ def make_image_temp( Parameters ---------- amr - ro ramses output object + ro pymses.RamsesOutput object center 3D array for coordinates center num output number @@ -317,7 +319,6 @@ def make_image_temp( ], origin="lower", ) - P.locator_params(axis="x", nbins=ntick) P.locator_params(axis="y", nbins=ntick) @@ -374,16 +375,13 @@ def disk_prop( directory_out = path_in # Check if the output file exists, and exit if it is the case - if ( - not force - and len(glob.glob(directory_out + "/prop_disk_" + str(num).zfill(5) + ".save")) - != 0 - ): + name_save = directory_out + "/prop_disk_" + str(num).zfill(5) + ".save" + if not force and len(glob.glob(name_save)) != 0: return # Compute the bins array lrad = np.log10(rad_ext) - rad = np.logspace(lrad - 2.0, lrad, num=nbin) + rad = np.logspace(lrad - 2.0, lrad, num=nb_bin) # Get Ramses data ro = pymses.RamsesOutput(path_in, num) @@ -439,7 +437,7 @@ def disk_prop( height = cs * np.sqrt(rc ** 3 / (G * mass_star)) mask_pos = np.abs(pos[:, 2]) < height # condition on position mask_dens = cells["rho"] > 1.0e6 # condition on density - mask = mask_pos or mask_dens + mask = mask_pos | mask_dens print("Number of selected cells ", np.sum(mask)) pos_disk = pos[mask] @@ -455,101 +453,72 @@ def disk_prop( # TODO Check what do that does nzoom = 9 eps = 0.5 ** nzoom - map_coldens, map_w13, xedges, yedges = me.make_hierarch_map( - pos_disk_x, - pos_disk_y, - pos_disk_z, - dx_disk, - rho_disk, - rho_disk, - eps, - center=[0.0, 0.0, 0.0], - make_image=do_plot, - path_out=directory_out, - tag="xy_" + str(num).zfill(5), - ) + # map_coldens , map_w13, xedges, yedges = me.make_hierarch_map(pos_disk_x,pos_disk_y,pos_disk_z,dx_disk,rho_disk,rho_disk,eps,center=[0.,0.,0.],make_image=do_plot,path_out=directory_out,tag='xy_'+ str(num).zfill(5)) - map_coldens, map_w13, xedges, yedges = me.make_hierarch_map( - pos_disk_z, - pos_disk_x, - pos_disk_y, - dx_disk, - rho_disk, - rho_disk, - eps, - center=[0.0, 0.0, 0.0], - make_image=do_plot, - path_out=directory_out, - tag="xz_" + str(num).zfill(5), - ) + # map_coldens , map_w13, xedges, yedges = me.make_hierarch_map(pos_disk_z,pos_disk_x,pos_disk_y,dx_disk,rho_disk,rho_disk,eps,center=[0.,0.,0.],make_image=do_plot,path_out=directory_out,tag='xz_'+ str(num).zfill(5)) # Initialize binned quantities norm_rad = lbox * scale_l / AU # radius in AU - rdisk_AU = rad_disk * norm_rad + rdisk_AU = rc_disk * norm_rad - cs_rad = np.zeros(nbin - 1) - temp_rad = np.zeros(nbin - 1) - press_rad = np.zeros(nbin - 1) - rho_rad = np.zeros(nbin - 1) - coldens_rad = np.zeros(nbin - 1) - v_az_rad = np.zeros(nbin - 1) - v_rad_rad = np.zeros(nbin - 1) + cs_rad = np.zeros(nb_bin - 1) + temp_rad = np.zeros(nb_bin - 1) + press_rad = np.zeros(nb_bin - 1) + rho_rad = np.zeros(nb_bin - 1) + coldens_rad = np.zeros(nb_bin - 1) + v_az_rad = np.zeros(nb_bin - 1) + v_rad_rad = np.zeros(nb_bin - 1) + alpha_rey_rad = np.zeros(nb_bin - 1) - for i in range(nbin - 1): - mask_bin = (rdisk_AU > rad[i]) or (rdisk_AU < rad[i + 1]) + for i in range(nb_bin - 1): + mask_bin = (rdisk_AU > rad[i]) & (rdisk_AU < rad[i + 1]) - press_rad[i] = np.sum(press_disk[mask] * dvol_disk[mask]) / np.sum( - dvol_disk[mask] + print( + "Bin {} cells between {} and {} AU".format( + np.sum(mask_bin), rad[i], rad[i + 1] + ) + ) + + press_rad[i] = np.sum(press_disk[mask_bin] * dvol_disk[mask_bin]) / np.sum( + dvol_disk[mask_bin] + ) + rho_rad[i] = np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) / np.sum( + dvol_disk[mask_bin] ) - rho_rad[i] = np.sum(rho_disk[mask] * dvol_disk[mask]) / np.sum(dvol_disk[mask]) temp_rad[i] = press_rad[i] / rho_rad[i] - # TODO vérifier unités + # TODO verifier unites # Surface of a bin : S = dr * 2 * pi * r with # dr = rad[i + 1] - rad[i] and r = (rad[i + 1] + rad[i]) / 2. coldens_rad[i] = ( - np.sum(rho_disk[mask] * dvol_disk[mask]) + np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) * (lbox * pc) ** 3 / ((rad[i + 1] - rad[i]) * (rad[i + 1] + rad[i]) * np.pi * AU ** 2) ) v_az_rad[i] = np.sum( - v_az_disk[mask] * rho_disk[mask] * dvol_disk[mask] - ) / np.sum(rho_disk[mask] * dvol_disk[mask]) + v_az_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin] + ) / np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) v_rad_rad[i] = np.sum( - v_rad_disk[mask] * rho_disk[mask] * dvol_disk[mask] - ) / np.sum(rho_disk[mask] * dvol_disk[mask]) + v_rad_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin] + ) / np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) alpha_rey_rad[i] = ( ( np.sum( - v_az_disk[mask] - * v_rad_disk[mask] - * rho_disk[mask] - * dvol_disk[mask] + v_az_disk[mask_bin] + * v_rad_disk[mask_bin] + * rho_disk[mask_bin] + * dvol_disk[mask_bin] ) - / np.sum(dvol_disk[mask] * press_disk[mask]) + / np.sum(dvol_disk[mask_bin] * press_disk[mask_bin]) - v_az_rad[i] * v_rad_rad[i] * rho_rad[i] / press_rad[i] ) * v_az_rad[i] / abs(v_az_rad[i]) ) - prop_disk = { - "time": time, - "rad_AU": rad[0 : nbin - 1], - "center": pos_mass, - "alpha_rey": alpha_rey_rad, - "v_rad": v_rad_rad, - "v_az": v_az_rad, - "coldens": coldens_rad, - "rho": rho_rad, - "press": press_rad, - "temp": temp_rad, - "cs": cs_rad, - } - # Convert to good units (TODO check) cs_rad = np.sqrt(temp_rad) * scale_v / km_s temp_rad = temp_rad * scale_T2 @@ -558,8 +527,21 @@ def disk_prop( v_az_rad = v_az_rad * scale_v / km_s v_rad_rad = v_rad_rad * scale_v / km_s + prop_disk = { + "time": time, + "rad_AU": rad[0 : nb_bin - 1], + "center": pos_star, + "alpha_rey": alpha_rey_rad, + "v_rad": v_rad_rad, + "v_az": v_az_rad, + "coldens": coldens_rad, + "rho": rho_rad, + "press": press_rad, + "temp": temp_rad, + "cs": cs_rad, + } + # store the results - name_save = directory_out + "/prop_disk_" + str(num).zfill(5) + ".save" f = open(name_save, "w") pickle.dump(prop_disk, f) f.close() @@ -582,16 +564,15 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): # Check if the properties file exists if len(glob.glob(name_save)) == 0: - print("no pickle file for disk properties. Run single_z_disc_prop") - return + throw("no pickle file for disk properties. Run single_z_disk_prop") f = open(name_save, "r") prop_disk = pickle.load(f) f.close() P.clf() P.plot( - np.log10(prop_disc["rad_AU"]), - np.log10(prop_disc["rho"]), + np.log10(prop_disk["rad_AU"]), + np.log10(prop_disk["rho"]), color="k", linewidth=2, ) @@ -599,31 +580,40 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): P.xlabel("disk radius (AU)") 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.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + ".pdf") + P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + ".jpeg") P.clf() P.plot( - np.log10(prop_disc["rad_AU"]), - np.log10(prop_disc["temp"]), + np.log10(prop_disk["rad_AU"]), + np.log10(prop_disk["temp"]), color="k", linewidth=2, ) P.ylabel(r"$\log(T) \, (K)$") - P.xlabel("disc radius (AU)") + P.xlabel("disk radius (AU)") 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.savefig(path + "/T_disk_r_" + str(num).zfill(5) + ".pdf") + P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + ".jpeg") P.clf() P.xscale("log") P.yscale("symlog", linthreshy=0.01) - P.plot((prop_disc["rad_AU"]), ((prop_disc["v_rad"])), color="k", linewidth=2) - P.plot((prop_disc["rad_AU"]), (abs(prop_disc["v_az"])), color="r", linewidth=2) - P.plot((prop_disc["rad_AU"]), ((prop_disc["cs"])), color="c", linewidth=2) + P.plot((prop_disk["rad_AU"]), ((prop_disk["v_rad"])), color="k", 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.legend((r"$v_r$", r"$v_\phi$", r"$c_s$"), loc="upper right") + + P.ylabel(r"$V \, (km s^{-1})$") + P.xlabel("disk radius (AU)") + + 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.legend((r"$v_r$", r"$v_\phi$", r"$c_s$"), loc="upper right") @@ -636,28 +626,32 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""): P.clf() P.plot( - np.log10(prop_disc["rad_AU"]), - np.log10(prop_disc["coldens"]), + np.log10(prop_disk["rad_AU"]), + np.log10(prop_disk["coldens"]), color="k", linewidth=2, ) P.ylabel(r"$\log(N) \, (cm^{-2})$") - P.xlabel("disc radius (AU)") + P.xlabel("disk radius (AU)") 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.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + ".pdf") + P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + ".jpeg") P.clf() P.xscale("log") P.yscale("symlog", linthreshy=0.001) + P.plot(prop_disk["rad_AU"], prop_disk["alpha_rey"], color="b", linewidth=2) + P.plot(prop_disc["rad_AU"], prop_disc["alpha_rey"], color="b", linewidth=2) - P.legend((r"$\alpha _{Rey}$"), loc="upper right") - P.ylabel(r"$\alpha}$") - P.xlabel("disc radius (AU)") + P.xlabel("disk radius (AU)") + + if pdf: + P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + ".pdf") + P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + ".jpeg") if pdf: P.savefig(path + "alpha_disk_r_" + str(num).zfill(5) + ".pdf") diff --git a/pipeline_disk.py b/pipeline_disk.py new file mode 100644 index 0000000..eedc592 --- /dev/null +++ b/pipeline_disk.py @@ -0,0 +1,58 @@ +# coding: utf-8 +import sys +import numpy as np +import module_extract as me +import extract_disk as d +import pymses +import os +import argparse +import disk_postprocess as dp + + +storage_in = "/drf/projets/alfven-data/" +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("-s", "--step", help="step between two output", type=int, default=1) +args = parser.parse_args() + + +user = "nbrucy" +folder = "simus" +project = "disk" +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 + path_out = storage_out + path_suffix + + if not os.path.exists(path_out): + os.makedirs(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, + ) + dp.disk_prop(path_in, i, path_out=path_out, rad_ext=50000) + dp.plot_disk_prop(path_out, i, tag=run + "_")