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

Add pdf for coldens

Browse Source
This commit is contained in:
Noe Brucy
2019-05-13 16:03:46 +02:00
parent 7b3793cff6
commit db17181fea
2 changed files with 373 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files
disk_postprocess.py
+325 -57
View File
@@ -1,9 +1,9 @@
# coding: utf-8
import sys
import numpy as np
import os
import pymses
import pymses
import numpy as np
import matplotlib as mpl
if os.environ.get("DISPLAY", "") == "":
@@ -12,7 +12,13 @@ if os.environ.get("DISPLAY", "") == "":
import pylab as P
import glob as glob
import pickle as pickle
try:
import cPickle as pickle
except:
print("cPickle not found, using pickle")
import pickle
import tables
from pymses.sources.ramses import output
from pymses.analysis import Camera, raytracing, slicing, splatting
@@ -31,12 +37,12 @@ def make_image_disk(
num,
path_out=None,
order="<",
save_data=True,
force=False,
tag="",
vel_red=20,
map_size=512,
put_title=True,
cpuamr=False,
cpu=False,
level=False,
pos_star=np.array([1.0, 1.0, 1.0]),
@@ -51,7 +57,7 @@ def make_image_disk(
path path of the Ramses output
num Ramses output number
path_out path of the pipeline outputb
order '<' or '>' TODO
order '<' or '>' order used by pymses for reading ramses 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
@@ -64,7 +70,7 @@ def make_image_disk(
rad = 0.5
center = [0.5, 0.5, 0.5]
make_image_aux(
return make_image_aux(
amr,
ro,
center,
@@ -73,10 +79,10 @@ def make_image_disk(
path,
force=force,
path_out=path_out,
save_data=save_data,
map_size=map_size,
vel_red=vel_red,
tag=tag,
cpuamr=cpuamr,
cpu=cpu,
level=level,
put_title=put_title,
@@ -95,10 +101,10 @@ def make_image_aux(
path,
force=False,
path_out=None,
save_data=True,
map_size=512,
vel_red=20,
tag="",
cpuamr=False,
cpu=False,
level=False,
pos_star=np.array([1.0, 1.0, 1.0]),
@@ -121,39 +127,39 @@ def make_image_aux(
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
"""
cpu = cpu or cpuamr
level = level or cpuamr
lbox = ro.info["boxlen"] # boxlen in codeunits
lbox_units = lbox
G = 1.0 # Gravitational constant
# Plotting parameters
ntick = 6
title_ax = {"x": "x (code)", "y": "y (code)", "z": "z (code)"}
im_extent = [
(-radius + center[0]) * lbox_units,
(radius + center[0]) * lbox_units,
(-radius + center[1]) * lbox_units,
(radius + center[1]) * lbox_units,
(-radius + center[0]) * lbox,
(radius + center[0]) * lbox,
(-radius + center[1]) * lbox,
(radius + center[1]) * lbox,
]
time = ro.info["time"] # time in codeunits
title = "t=" + str(time)[0:5] + " (code)"
# Determining outpout directory
if path_out is not None:
directory = path_out
else:
directory = path
# Checking for existing files
name = directory + "/coldens_z" + "_" + tag + "_" + format(num, "05") + out_ext
if len(glob.glob(name)) == 1 and not force:
return
# Prepare saving data
if save_data:
maps_disk = {"time": time, "im_extent": im_extent}
rho_op = ScalarOperator(lambda dset: dset["rho"], ro.info["unit_density"])
rt = None
@@ -189,15 +195,19 @@ def make_image_aux(
# Levels
if level:
level_op = MaxLevelOperator()
amr.set_read_levelmax(20)
level_op = MaxLevelOperator()
rt_level = raytracing.RayTracer(amr, ro.info, level_op)
datamap = rt_level.process(cam, surf_qty=True)
map_level = datamap.map.T
# if save_data:
# maps_disk['levels_' + ax_los] = map_level
levels_ar = np.arange(ro.info["levelmin"], ro.info["levelmax"] + 1)
# Computing linewidths
lw = np.ones(levels_ar.size) * 2
lvl_th = 8
lvl_th = 8 # Level threeshold for reducing linewidths
lw[levels_ar >= lvl_th] = lw[levels_ar >= lvl_th] ** (
lvl_th - levels_ar[levels_ar >= lvl_th]
)
@@ -226,6 +236,9 @@ def make_image_aux(
dmap_col = datamap.map.T * lbox
map_col = np.log10(dmap_col)
if save_data:
maps_disk["coldens_" + ax_los] = dmap_col
im = P.imshow(map_col, extent=im_extent, origin="lower")
P.locator_params(axis=ax_h, nbins=ntick)
@@ -249,9 +262,9 @@ def make_image_aux(
P.close()
# 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
datamap_rho = slicing.SliceMap(amr, cam, rho_op, z=0.0)
dmap_rho = (datamap_rho.map).T
map_rho = np.log10(dmap_rho)
vh_op = ScalarOperator(
lambda dset: dset["vel"][:, ax_nb[ax_h]], ro.info["unit_velocity"]
@@ -260,7 +273,6 @@ def make_image_aux(
map_vh_red = dmap_vh.map[
::vel_red, ::vel_red
] # take only a subset of velocities
map_vh_red = map_vh_red.T
vv_op = ScalarOperator(
@@ -270,6 +282,11 @@ def make_image_aux(
map_vv_red = dmap_vv.map[::vel_red, ::vel_red]
map_vv_red = map_vv_red.T
# if save_data:
# maps_disk['rho_' + ax_los] = dmap_rho
# maps_disk['v' + ax_h + '_' + ax_los] = (dmap_vh.map).T
# maps_disk['v' + ax_v + '_' + ax_los] = (dmap_vv.map).T
im = P.imshow(map_rho, extent=im_extent, origin="lower")
P.locator_params(axis=ax_h, nbins=ntick)
P.locator_params(axis=ax_v, nbins=ntick)
@@ -277,12 +294,13 @@ def make_image_aux(
nv = map_vv_red.shape[1]
vec_h = (
np.arange(nh) * 2.0 / nh * radius - radius + center[0] + radius / nh
) * lbox_units
) * lbox
vec_v = (
np.arange(nv) * 2.0 / nv * radius - radius + center[1] + radius / nv
) * lbox_units
) * lbox
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,
@@ -293,7 +311,6 @@ def make_image_aux(
labelpos="E",
coordinates="figure",
)
if put_title:
P.title(title)
P.xlabel(title_ax[ax_h])
@@ -310,11 +327,14 @@ def make_image_aux(
P.close()
P_op = ScalarOperator(lambda dset: dset["P"], ro.info["unit_pressure"])
dmap_P = slicing.SliceMap(amr, cam, P_op, z=0.0)
dmap_P = (slicing.SliceMap(amr, cam, P_op, z=0.0)).map.T
dmap_T = dmap_P.map.T / dmap_rho.map.T
dmap_T = dmap_P / dmap_rho
map_T = np.log10(dmap_T)
# if save_data:
# maps_disk['T_' + ax_los] = dmap_T
im = P.imshow(map_T, extent=im_extent, origin="lower")
P.locator_params(axis="x", nbins=ntick)
@@ -372,8 +392,10 @@ def make_image_aux(
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 = dmap_Q
map_Q = (lbox * dmap_cs.map.T) * dmap_omega.map.T / (np.pi * G * dmap_col)
# if save_data:
# maps_disk['Q_' + ax_los] = map_Q
im = P.imshow(
map_Q, extent=im_extent, origin="lower", norm=mpl.colors.LogNorm()
@@ -406,6 +428,9 @@ def make_image_aux(
datamap = rt_cpu.process(cam, surf_qty=True)
map_cpu = datamap.map.T
# if save_data:
# maps_disk['cpu_' + ax_los] = map_cpu
im = P.imshow(map_cpu, extent=im_extent, origin="lower")
P.locator_params(axis="x", nbins=ntick)
P.locator_params(axis="y", nbins=ntick)
@@ -427,6 +452,15 @@ def make_image_aux(
P.savefig(name_im)
P.close()
if save_data:
name_save = (
directory + "/maps_disk" + "_" + tag + "_" + format(num, "05") + ".save"
)
f = open(name_save, "w")
pickle.dump(maps_disk, f)
f.close()
return maps_disk
def disk_prop(
path_in,
@@ -437,6 +471,7 @@ def disk_prop(
rad_ext=1.0,
mass_star=1.0,
pos_star=np.array([1.0, 1.0, 1.0]),
binning="log",
):
"""
Compute properties of a disk in the plane (0,x,y)
@@ -470,11 +505,16 @@ def disk_prop(
if not force and len(glob.glob(name_save)) != 0:
return
nb_bin_hist = nb_bin
# Compute the bins array
if binning == "log":
lrad = np.log10(rad_ext)
rad = np.logspace(lrad - 2.0, lrad, num=nb_bin)
elif binning == "lin":
rad = np.linspace(0.0, rad_ext, num=nb_bin)
else:
raise ValueError(
"Incorrect binning specification, binnning should be 'lin' or 'log'"
)
# Get Ramses data
ro = pymses.RamsesOutput(path_in, num)
@@ -553,10 +593,6 @@ def disk_prop(
surf_rad = np.zeros(nb_bin - 1) # Surface of a bin
mass_rad = np.zeros(nb_bin - 1) # Mass of a bin
# Density fluctuations
hist_drho = np.zeros(nb_bin_hist)
hist_edges = np.zeros(nb_bin_hist + 1)
for i in range(nb_bin - 1):
mask_bin = (rc_disk > rad[i]) & (rc_disk < rad[i + 1])
@@ -620,20 +656,12 @@ def disk_prop(
/ mass_rad[i]
)
# Histogramm : density fluctuaction distribution function
drho = np.log(rho_disk[mask_bin] / rho_rad[i])
hist, hist_edges = P.histogram(
drho, bins=nb_bin_hist, weights=dvol_disk[mask_bin]
)
hist_drho = hist_drho + hist
# Derived quantities
cs_rad = np.sqrt(temp_rad)
Q_kepl_rad = cs_rad * v_az_rad / (np.pi * G * coldens_rad * rad[0 : nb_bin - 1])
# Means
mask_mean = (0.1 < rad[0 : nb_bin - 1]) & (rad[0 : nb_bin - 1] < 0.2)
print(rad[0 : nb_bin - 1][mask_mean])
mass_mean = np.sum(mass_rad[mask_mean])
alpha_rey_mean = np.sum(alpha_rey_rad[mask_mean] * mass_rad[mask_mean]) / mass_mean
alpha_grav_mean = (
@@ -641,7 +669,6 @@ def disk_prop(
)
Q_mean = np.sum(Q_kepl_rad[mask_mean] * mass_rad[mask_mean]) / mass_mean
Q_min = np.nanmin(Q_kepl_rad)
print("alphas, Q ", alpha_rey_mean, alpha_grav_mean, Q_mean)
# store the results
prop_disk = {
@@ -660,8 +687,6 @@ def disk_prop(
"coldens": coldens_rad,
"rho": rho_rad,
"press": press_rad,
"hist_drho": hist_drho,
"hist_edges": hist_edges,
"temp": temp_rad,
"cs": cs_rad,
"Q_kepl": Q_kepl_rad,
@@ -837,21 +862,105 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False):
P.savefig(path + "/H_r_" + str(num).zfill(5) + out_ext)
P.close()
# Density fluctuation histogram
def disk_pdf(
path,
num,
maps_disk,
pos_star=[1.0, 1.0],
force=False,
interactive=False,
nb_bin_hist=50,
tag="",
rad_min=0.1,
):
# Load property file
name_prop = path + "/prop_disk_" + str(num).zfill(5) + ".save"
# Check if the properties file exists
if len(glob.glob(name_prop)) == 0:
raise IOError("no pickle file for disk properties. Run disk_prop() first")
f = open(name_prop, "r")
prop_disk = pickle.load(f)
f.close()
# Load maps file
print("load maps file")
name_maps = path + "/maps_disk" + "_" + tag + "_" + format(num, "05") + ".save"
# Check if the maps file exists
if maps_disk is None:
if len(glob.glob(name_maps)) == 0:
raise IOError("no pickle file for disk maps. Run make_image_disk() first")
f = open(name_maps, "r")
maps_disk = pickle.load(f)
f.close()
print("maps file loaded")
time = prop_disk["time"]
title = "t=" + str(time)[0:5] + " (code)"
coldens_map = maps_disk["coldens_z"]
im_extent = maps_disk["im_extent"]
rad_bins = prop_disk["rad"]
coldens_rad = prop_disk["coldens"]
x = np.linspace(im_extent[0], im_extent[1], coldens_map.shape[0])
y = np.linspace(im_extent[2], im_extent[3], coldens_map.shape[1])
xx, yy = np.meshgrid(x, y)
rr = np.sqrt((xx - pos_star[0]) ** 2 + (yy - pos_star[1]) ** 2)
# Find appropriate bin for each coordinate set
bins = np.zeros(rr.shape, dtype=int)
for r in rad_bins[1:]:
bins = bins + (rr >= r).astype(int)
coldens_mean = coldens_rad[bins.flatten()]
coldens_mean_map = np.reshape(coldens_mean, coldens_map.shape)
# Kill fluctuations for r < rad_min
coldens_map[rr < rad_min] = coldens_mean_map[rr < rad_min]
# Histogramm : density fluctuation distribution function
dcoldens = np.log10(coldens_map.flatten() / coldens_mean)
P.grid()
P.xlabel(r"$\log(\frac{\rho}{\bar{\rho}})$")
P.ylabel(r"fraction of total volume")
P.yscale("log")
P.xlabel(r"$\log(N / \bar{N})$")
P.ylabel("Probability density")
P.title(title)
hist = prop_disk["hist_drho"]
egdes = prop_disk["hist_edges"]
widths = egdes[1:] - egdes[:-1]
centers = egdes[:-1] + widths / 2.0
P.bar(centers, hist, width=widths)
P.hist(dcoldens, nb_bins_hist, range=(-1, 2))
if interactive:
pass
else:
P.savefig(path + "/drho_hist_" + str(num).zfill(5) + out_ext)
P.savefig(path + "/dcol_hist_" + str(num).zfill(5) + out_ext)
P.close()
# Fluctuation map
dcoldens_map = np.reshape(dcoldens, coldens_map.shape)
# cont = P.contour(coldens_mean_map,
# extent=im_extent,
# origin='lower',
# colors='k',
# linewidths=0.1)
im = P.imshow(dcoldens_map, extent=im_extent, origin="lower", cmap="viridis")
P.title(title)
P.xlabel(r"$x$")
P.ylabel(r"$y$")
cbar = P.colorbar(im)
cbar.set_label(r"$log(N/\bar{N})$ (code)")
name = path + "/dcoldensmap_" + format(num, "05")
name_im = name + out_ext
if interactive:
P.figure()
else:
P.savefig(name_im)
P.close()
@@ -1026,3 +1135,162 @@ def evolution(path, nums, force=False, interactive=False):
else:
P.savefig(path + "/mass_time" + out_ext)
P.close()
def make_clump_hop(
path_in,
num,
name,
thres_dens,
thres_level,
pos_zoom,
size_zoom,
path_out=None,
path_hop="",
force=False,
gcomp=True,
):
"""
This routine use the HOP algorithm to extract clumps defined from their peaks
as an output it provides a list of cell position ordered by the group to which they belong
Parameters
----------
path_in is the path where the data tobe read are located
path_out is the path of teh directory where resulting files must be written
num output number
name a string which is used to write the names of the various files
thres_dens density threshold above which cells are considered
thres_level level threshold above which cells are considered
pos_zoom the center of the zoom coordinates
size_zoom the 3 zoom extension (x, y and z)
"""
if path_out is not None:
directory_out = path_out
else:
directory_out = path_in
name_txt = name + ".txt"
# check whether hop entry files have been created (not test is done on .txt only
if len(glob.glob(name_txt)) == 0 or force:
ro = pymses.RamsesOutput(path_in, num)
amr = ro.amr_source(["rho"], grav_compat=gcomp) # density only is used
center = np.asarray(pos_zoom)
radius = size_zoom
min_coords = np.zeros(3)
max_coords = np.zeros(3)
min_coords[:] = center[:] - radius / 2.0
max_coords[:] = center[:] + radius / 2.0
region = Box((min_coords, max_coords))
# region = Sphere(center,radius)
filt_amr = RegionFilter(region, amr)
cell_source = CellsToPoints(
filt_amr,
)
# selection of the cells of interest
def cell_selec_func(dset):
mask1 = dset["rho"] >= thres_dens
dx = dset.get_sizes()
mask2 = dx <= 0.5 ** thres_level
return mask1 * mask2
# begin cell_selec
cells_selec = PointFunctionFilter(cell_selec_func, cell_source).flatten()
dx = cells_selec.get_sizes()
ncells = cells_selec.npoints
# fill the matrice with ID, x,y,z and masses of particles
val_mat = np.zeros((ncells, 5))
val_mat[:, 0] = np.arange(ncells)
val_mat[:, 1:4] = cells_selec.points[:, 0:3]
val_mat[:, 4] = cells_selec["rho"] * (dx ** 3)
# write name.txt
head = str(ncells)
np.savetxt(
name_txt,
val_mat,
fmt="%10d %.10e %.10e %.10e %.10e",
header=head,
delimiter=" ",
comments=" ",
)
# end of creation name.txt
# creation name.den
f = open(name + ".den", "wb")
f.write(pack("I", ncells))
cells_selec["rho"].astype("f").tofile(f)
f.close()
print(name + ".den created")
# end of creation name.den
# HOP Algorithm
print("creation of .hop and .gbound du to hop")
fname = path_hop + name + ".txt"
print("look for hop in ", fname)
h = HOP(fname, path_hop)
h.process_hop()
print("hop grouping is finished")
# end of HOP Algorithm
idpart = val_mat[:, 0]
X = val_mat[:, 1]
Y = val_mat[:, 2]
Z = val_mat[:, 3]
mass = val_mat[:, 4]
# read the gbound file to get list of particle numbers within groups
f = open(name + ".gbound", "r")
aline = f.readline()
ngroups = int(aline)
npart_v = np.zeros(ngroups, dtype=int)
for i in range(10):
aline = f.readline()
for i in range(ngroups):
aline = f.readline()
vec = aline.split()
igroup = int(vec[0])
npart_v[igroup] = int(vec[1])
f.close()
# get the igroup array
group_ids = h.get_group_ids()
# sort it and apply the sorting to the coordinates
# this means that the particules of group 1 are written first then of group 2 etc...
ind_sort = np.argsort(group_ids)
xx_v = X[ind_sort]
yy_v = Y[ind_sort]
zz_v = Z[ind_sort]
vect_id_group = group_ids[ind_sort] # not so useful
# write the sorted cells
name_save_clump = directory_out + name + ".save"
np.savez(
name_save_clump,
ngroups=ngroups,
npart_v=npart_v,
xx_v=xx_v,
yy_v=yy_v,
zz_v=zz_v,
vect_id_group=vect_id_group,
num=num,
name=name,
thres_dens=thres_dens,
thres_level=thres_level,
pos_zoom=pos_zoom,
size_zoom=size_zoom,
)
return name_save_clump
pipeline_disk.py
+46 -8
View File
@@ -57,17 +57,36 @@ parser.add_argument(
help="plot evolution of quantities over time",
action="store_true",
)
parser.add_argument(
"--pdf", help="plot pdf of fluctuations of column density", action="store_true"
)
parser.add_argument(
"--fft", help="use quick and dirty fft rendering", action="store_true"
)
parser.add_argument("--level", help="plot levels", action="store_true")
parser.add_argument("--cpu", help="plot cpu", action="store_true")
parser.add_argument(
"-ms",
"--mapsize",
help="size of the maps created in he map mode (in pixel)",
type=int,
default=1024,
)
parser.add_argument(
"--nb_bin", help="Number of bins for azimuthal averages", type=int, default=50
)
parser.add_argument(
"--binning",
help="Kind of binning (logarithmic or linear)",
choices=["log", "lin"],
default="log",
)
parser.add_argument(
"--rad_ext", help="Value of the highest bin", type=float, default=1.0
)
parser.add_argument(
"-x", help="x position of the central point", type=float, default=1.0
)
@@ -125,13 +144,15 @@ for run in runs:
while not success:
try:
maps_disk = None
if args.maps:
dp.make_image_disk(
print("[{}, {}] computing maps".format(run, i))
maps_disk = dp.make_image_disk(
path_in,
i,
path_out=path_out,
tag=run,
map_size=1024,
map_size=args.mapsize,
force=args.force_redo,
level=args.level,
cpu=args.cpu,
@@ -139,17 +160,22 @@ for run in runs:
fft=args.fft,
interactive=args.interactive,
)
if args.disk or args.compare or args.evolution:
print("[{}, {}] maps computed".format(run, i))
if args.disk:
print("[{}, {}] computing disk props".format(run, i))
dp.disk_prop(
path_in,
i,
path_out=path_out,
rad_ext=1,
nb_bin=args.nb_bin,
binning=args.binning,
rad_ext=args.rad_ext,
force=args.force_redo,
pos_star=np.array([args.x, args.y, args.z]),
)
print("[{}, {}] disk_props computed".format(run, i))
if args.disk:
print("[{}, {}] plotting disk props".format(run, i))
dp.plot_disk_prop(
path_out,
i,
@@ -157,14 +183,26 @@ for run in runs:
force=args.force_redo,
interactive=args.interactive,
)
print("[{}, {}] disk props plotted".format(run, i))
if args.pdf:
print("[{}, {}] computing pdf".format(run, i))
dp.disk_pdf(
path_out,
i,
maps_disk,
pos_star=np.array([args.x, args.y, args.z]),
force=args.force_redo,
tag=run,
interactive=args.interactive,
)
print("[{}, {}] pdf computed".format(run, i))
success = True
except ValueError:
except ValueError as e:
print(e)
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(
"Unable to proceed for run {} output {}. Trying again in {} s ({} tries remaining)".format(
run, i, args.waiting_time, args.allowed_failures - failures
)
)