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

adapt postprocess

Browse Source
This commit is contained in:
Noe Brucy
2019-04-05 09:46:24 +02:00
parent b852cfec80
commit dca8ab5d41
2 changed files with 149 additions and 97 deletions
Download Patch File Download Diff File Expand all files Collapse all files
disk_postprocess.py
+87 -93
View File
@@ -1,11 +1,13 @@
# coding: utf-8
import sys import sys
import numpy as np import numpy as np
import os import os
import pymses import pymses
import pylab as P import pylab as P
import glob as glob import glob as glob
import pandas as pd
import pickle as pickle import pickle as pickle
import module_extract as me
from pymses.filters import CellsToPoints
from pymses.sources.ramses import output from pymses.sources.ramses import output
@@ -38,7 +40,7 @@ def make_image_temp(
Parameters Parameters
---------- ----------
amr amr
ro ramses output object ro pymses.RamsesOutput object
center 3D array for coordinates center center 3D array for coordinates center
num output number num output number
@@ -317,7 +319,6 @@ def make_image_temp(
], ],
origin="lower", origin="lower",
) )
P.locator_params(axis="x", nbins=ntick) P.locator_params(axis="x", nbins=ntick)
P.locator_params(axis="y", nbins=ntick) P.locator_params(axis="y", nbins=ntick)
@@ -374,16 +375,13 @@ def disk_prop(
directory_out = path_in directory_out = path_in
# Check if the output file exists, and exit if it is the case # Check if the output file exists, and exit if it is the case
if ( name_save = directory_out + "/prop_disk_" + str(num).zfill(5) + ".save"
not force if not force and len(glob.glob(name_save)) != 0:
and len(glob.glob(directory_out + "/prop_disk_" + str(num).zfill(5) + ".save"))
!= 0
):
return return
# Compute the bins array # Compute the bins array
lrad = np.log10(rad_ext) 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 # Get Ramses data
ro = pymses.RamsesOutput(path_in, num) ro = pymses.RamsesOutput(path_in, num)
@@ -439,7 +437,7 @@ def disk_prop(
height = cs * np.sqrt(rc ** 3 / (G * mass_star)) height = cs * np.sqrt(rc ** 3 / (G * mass_star))
mask_pos = np.abs(pos[:, 2]) < height # condition on position mask_pos = np.abs(pos[:, 2]) < height # condition on position
mask_dens = cells["rho"] > 1.0e6 # condition on density 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)) print("Number of selected cells ", np.sum(mask))
pos_disk = pos[mask] pos_disk = pos[mask]
@@ -455,91 +453,84 @@ def disk_prop(
# TODO Check what do that does # TODO Check what do that does
nzoom = 9 nzoom = 9
eps = 0.5 ** nzoom eps = 0.5 ** nzoom
map_coldens, map_w13, xedges, yedges = me.make_hierarch_map( # 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))
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( # 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))
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),
)
# Initialize binned quantities # Initialize binned quantities
norm_rad = lbox * scale_l / AU # radius in AU 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) cs_rad = np.zeros(nb_bin - 1)
temp_rad = np.zeros(nbin - 1) temp_rad = np.zeros(nb_bin - 1)
press_rad = np.zeros(nbin - 1) press_rad = np.zeros(nb_bin - 1)
rho_rad = np.zeros(nbin - 1) rho_rad = np.zeros(nb_bin - 1)
coldens_rad = np.zeros(nbin - 1) coldens_rad = np.zeros(nb_bin - 1)
v_az_rad = np.zeros(nbin - 1) v_az_rad = np.zeros(nb_bin - 1)
v_rad_rad = np.zeros(nbin - 1) v_rad_rad = np.zeros(nb_bin - 1)
alpha_rey_rad = np.zeros(nb_bin - 1)
for i in range(nbin - 1): for i in range(nb_bin - 1):
mask_bin = (rdisk_AU > rad[i]) or (rdisk_AU < rad[i + 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( print(
dvol_disk[mask] "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] 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 # Surface of a bin : S = dr * 2 * pi * r with
# dr = rad[i + 1] - rad[i] and r = (rad[i + 1] + rad[i]) / 2. # dr = rad[i + 1] - rad[i] and r = (rad[i + 1] + rad[i]) / 2.
coldens_rad[i] = ( coldens_rad[i] = (
np.sum(rho_disk[mask] * dvol_disk[mask]) np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin])
* (lbox * pc) ** 3 * (lbox * pc) ** 3
/ ((rad[i + 1] - rad[i]) * (rad[i + 1] + rad[i]) * np.pi * AU ** 2) / ((rad[i + 1] - rad[i]) * (rad[i + 1] + rad[i]) * np.pi * AU ** 2)
) )
v_az_rad[i] = np.sum( v_az_rad[i] = np.sum(
v_az_disk[mask] * rho_disk[mask] * dvol_disk[mask] v_az_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin]
) / np.sum(rho_disk[mask] * dvol_disk[mask]) ) / np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin])
v_rad_rad[i] = np.sum( v_rad_rad[i] = np.sum(
v_rad_disk[mask] * rho_disk[mask] * dvol_disk[mask] v_rad_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin]
) / np.sum(rho_disk[mask] * dvol_disk[mask]) ) / np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin])
alpha_rey_rad[i] = ( alpha_rey_rad[i] = (
( (
np.sum( np.sum(
v_az_disk[mask] v_az_disk[mask_bin]
* v_rad_disk[mask] * v_rad_disk[mask_bin]
* rho_disk[mask] * rho_disk[mask_bin]
* dvol_disk[mask] * 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] * v_rad_rad[i] * rho_rad[i] / press_rad[i]
) )
* v_az_rad[i] * v_az_rad[i]
/ abs(v_az_rad[i]) / 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_az_rad = v_az_rad * scale_v / km_s
v_rad_rad = v_rad_rad * scale_v / km_s
prop_disk = { prop_disk = {
"time": time, "time": time,
"rad_AU": rad[0 : nbin - 1], "rad_AU": rad[0 : nb_bin - 1],
"center": pos_mass, "center": pos_star,
"alpha_rey": alpha_rey_rad, "alpha_rey": alpha_rey_rad,
"v_rad": v_rad_rad, "v_rad": v_rad_rad,
"v_az": v_az_rad, "v_az": v_az_rad,
@@ -550,16 +541,7 @@ def disk_prop(
"cs": cs_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
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
# store the results # store the results
name_save = directory_out + "/prop_disk_" + str(num).zfill(5) + ".save"
f = open(name_save, "w") f = open(name_save, "w")
pickle.dump(prop_disk, f) pickle.dump(prop_disk, f)
f.close() f.close()
@@ -582,16 +564,15 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
# Check if the properties file exists # Check if the properties file exists
if len(glob.glob(name_save)) == 0: if len(glob.glob(name_save)) == 0:
print("no pickle file for disk properties. Run single_z_disc_prop") throw("no pickle file for disk properties. Run single_z_disk_prop")
return
f = open(name_save, "r") f = open(name_save, "r")
prop_disk = pickle.load(f) prop_disk = pickle.load(f)
f.close() f.close()
P.clf() P.clf()
P.plot( P.plot(
np.log10(prop_disc["rad_AU"]), np.log10(prop_disk["rad_AU"]),
np.log10(prop_disc["rho"]), np.log10(prop_disk["rho"]),
color="k", color="k",
linewidth=2, linewidth=2,
) )
@@ -599,31 +580,40 @@ def plot_disk_prop(path, num, force=False, pdf=False, tag=""):
P.xlabel("disk radius (AU)") P.xlabel("disk radius (AU)")
if pdf: if pdf:
P.savefig(path + "rho_disk_r_" + str(num).zfill(5) + ".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) + ".jpeg")
P.clf() P.clf()
P.plot( P.plot(
np.log10(prop_disc["rad_AU"]), np.log10(prop_disk["rad_AU"]),
np.log10(prop_disc["temp"]), np.log10(prop_disk["temp"]),
color="k", color="k",
linewidth=2, linewidth=2,
) )
P.ylabel(r"$\log(T) \, (K)$") P.ylabel(r"$\log(T) \, (K)$")
P.xlabel("disc radius (AU)") P.xlabel("disk radius (AU)")
if pdf: if pdf:
P.savefig(path + "T_disk_r_" + str(num).zfill(5) + ".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) + ".jpeg")
P.clf() P.clf()
P.xscale("log") P.xscale("log")
P.yscale("symlog", linthreshy=0.01) P.yscale("symlog", linthreshy=0.01)
P.plot((prop_disc["rad_AU"]), ((prop_disc["v_rad"])), color="k", linewidth=2) P.plot((prop_disk["rad_AU"]), ((prop_disk["v_rad"])), color="k", linewidth=2)
P.plot((prop_disc["rad_AU"]), (abs(prop_disc["v_az"])), color="r", linewidth=2) P.plot((prop_disk["rad_AU"]), (abs(prop_disk["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["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") 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.clf()
P.plot( P.plot(
np.log10(prop_disc["rad_AU"]), np.log10(prop_disk["rad_AU"]),
np.log10(prop_disc["coldens"]), np.log10(prop_disk["coldens"]),
color="k", color="k",
linewidth=2, linewidth=2,
) )
P.ylabel(r"$\log(N) \, (cm^{-2})$") P.ylabel(r"$\log(N) \, (cm^{-2})$")
P.xlabel("disc radius (AU)") P.xlabel("disk radius (AU)")
if pdf: if pdf:
P.savefig(path + "coldens_disk_r_" + str(num).zfill(5) + ".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) + ".jpeg")
P.clf() P.clf()
P.xscale("log") P.xscale("log")
P.yscale("symlog", linthreshy=0.001) 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.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.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: if pdf:
P.savefig(path + "alpha_disk_r_" + str(num).zfill(5) + ".pdf") P.savefig(path + "alpha_disk_r_" + str(num).zfill(5) + ".pdf")
pipeline_disk.py
+58
View File
@@ -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 + "_")