noe/pipeline
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

adapt disk analysis

Browse Source
This commit is contained in:
Noe Brucy
2019-04-08 10:54:31 +02:00
parent dca8ab5d41
commit 260dcfc7a5
2 changed files with 46 additions and 70 deletions
Download Patch File Download Diff File Expand all files Collapse all files
disk_postprocess.py
+36 -68
View File
@@ -9,7 +9,6 @@ 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
@@ -343,7 +342,7 @@ def disk_prop(
path_out=None,
force=False,
nb_bin=20,
rad_ext=100.0,
rad_ext=1.0,
mass_star=1.0,
pos_star=np.array([1.0, 1.0, 1.0]),
):
@@ -387,42 +386,21 @@ def disk_prop(
ro = pymses.RamsesOutput(path_in, num)
lbox = ro.info["boxlen"] # boxlen in codeunits (=>pc)
(
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,
) = me.normalisation(ro)
time = ro.info["time"] * scale_t / Myr
time = ro.info["time"] # * scale_t / Myr
# Get array of cell positions
amr = ro.amr_source(["rho", "vel", "Br", "Bl", "P"])
cell_source = CellsToPoints(amr)
cells = cell_source.flatten()
dx = cells.get_sizes()
pos = cells.points
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
vel = cells["vel"]
vel = cells["vel"] * lbox
# Get radial component of velocity
norm_pos = rc
norm_pos[rc == 0] = 1.0e-10 # Avoid division per 0
@@ -432,12 +410,12 @@ def disk_prop(
# Select cells that are actually in the disk, ie within the scale height
# TODO Check units
G = 6.8e-8
cs = np.sqrt(cells["P"] / cells["rho"]) * scale_v # sound velocity
G = 1.0 # G=6.8e-8
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
mask_dens = cells["rho"] > 1.0e6 # condition on density
mask = mask_pos | mask_dens
mask = mask_pos # & mask_dens
print("Number of selected cells ", np.sum(mask))
pos_disk = pos[mask]
@@ -449,33 +427,25 @@ def disk_prop(
dvol_disk = dx_disk ** 3
v_rad_disk = v_rad[mask]
v_az_disk = v_az[mask]
# 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.],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.],make_image=do_plot,path_out=directory_out,tag='xz_'+ str(num).zfill(5))
v_kepl = np.sqrt(mass_star * G / rc_disk)
# Initialize binned quantities
norm_rad = lbox * scale_l / AU # radius in AU
rdisk_AU = rc_disk * norm_rad
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_kepl_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(nb_bin - 1):
mask_bin = (rdisk_AU > rad[i]) & (rdisk_AU < rad[i + 1])
mask_bin = (rc_disk > rad[i]) & (rc_disk < rad[i + 1])
print(
"Bin {} cells between {} and {} AU".format(
np.sum(mask_bin), rad[i], rad[i + 1]
"Bin #{} : {} cells between {} and {}".format(
i, np.sum(mask_bin), rad[i], rad[i + 1]
)
)
@@ -492,8 +462,8 @@ def disk_prop(
# dr = rad[i + 1] - rad[i] and r = (rad[i + 1] + rad[i]) / 2.
coldens_rad[i] = (
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)
* (lbox) ** 3
/ ((rad[i + 1] - rad[i]) * (rad[i + 1] + rad[i]) * np.pi)
)
v_az_rad[i] = np.sum(
@@ -519,13 +489,18 @@ def disk_prop(
/ abs(v_az_rad[i])
)
# Convert to good units (TODO check)
cs_rad = np.sqrt(temp_rad) * scale_v / km_s
temp_rad = temp_rad * scale_T2
press_rad = press_rad * scale_v ** 2 * scale_d
v_kepl_rad[i] = np.sum(
v_kepl[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin]
) / np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin])
v_az_rad = v_az_rad * scale_v / km_s
v_rad_rad = v_rad_rad * scale_v / km_s
# Convert to good units (TODO check)
cs_rad = np.sqrt(temp_rad) # *scale_v / km_s
temp_rad = temp_rad # * scale_T2
press_rad = press_rad # * scale_v**2 * scale_d
v_az_rad = v_az_rad # * scale_v / km_s
v_rad_rad = v_rad_rad # * scale_v / km_s
v_kepl_rad = v_kepl_rad
prop_disk = {
"time": time,
@@ -534,6 +509,7 @@ def disk_prop(
"alpha_rey": alpha_rey_rad,
"v_rad": v_rad_rad,
"v_az": v_az_rad,
"v_kepl": v_kepl_rad,
"coldens": coldens_rad,
"rho": rho_rad,
"press": press_rad,
@@ -577,7 +553,7 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
linewidth=2,
)
P.ylabel(r"$\log(n) \, (cm^{-3})$")
P.xlabel("disk radius (AU)")
P.xlabel("disk radius")
if pdf:
P.savefig(path + "/rho_disk_r_" + str(num).zfill(5) + ".pdf")
@@ -591,7 +567,7 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
linewidth=2,
)
P.ylabel(r"$\log(T) \, (K)$")
P.xlabel("disk radius (AU)")
P.xlabel("disk radius")
if pdf:
P.savefig(path + "/T_disk_r_" + str(num).zfill(5) + ".pdf")
@@ -603,27 +579,19 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
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.legend((r"$v_r$", r"$v_\phi$", r"$c_s$"), loc="upper right")
P.legend((r"$v_r$", r"$v_{kepl}$", r"$v_\phi$", r"$c_s$"), loc="upper right")
P.ylabel(r"$V \, (km s^{-1})$")
P.xlabel("disk radius (AU)")
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.legend((r"$v_r$", r"$v_\phi$", r"$c_s$"), loc="upper right")
P.ylabel(r"$V \, (km s^{-1})$")
P.xlabel("disc 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.clf()
P.plot(
np.log10(prop_disk["rad_AU"]),
@@ -632,7 +600,7 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
linewidth=2,
)
P.ylabel(r"$\log(N) \, (cm^{-2})$")
P.xlabel("disk radius (AU)")
P.xlabel("disk radius ")
if pdf:
P.savefig(path + "/coldens_disk_r_" + str(num).zfill(5) + ".pdf")
@@ -642,12 +610,12 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
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_disk["rad_AU"], abs(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("disk radius (AU)")
P.ylabel(r"$\alpha$")
P.xlabel("disk radius ")
if pdf:
P.savefig(path + "/alpha_disk_r_" + str(num).zfill(5) + ".pdf")
pipeline_disk.py
+10 -2
View File
@@ -21,6 +21,10 @@ 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)
parser.add_argument(
"-d", "--disk", help="do specific disk radial analysis", action="store_true"
)
args = parser.parse_args()
@@ -54,5 +58,9 @@ for run in runs:
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 + "_")
# 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 + "_")