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

Corrected errors

Browse Source
This commit is contained in:
Noe Brucy
2019-06-14 14:14:02 +02:00
parent 84fcf81dfc
commit 3afc117e2a
2 changed files with 166 additions and 165 deletions
Download Patch File Download Diff File Expand all files Collapse all files
disk_postprocess.py
+148 -147
View File
@@ -21,6 +21,7 @@ except:
import tables
from scipy.stats import linregress
from pymses.sources.ramses import output
from pymses.sources.hop.file_formats import *
from pymses.analysis import Camera, raytracing, slicing, splatting
from pymses.filters import CellsToPoints
from pymses.analysis import ScalarOperator, FractionOperator, MaxLevelOperator
@@ -317,18 +318,15 @@ def plot_maps(
):
continue
map_disk = maps_disk[image + "_" + ax_los]
map_disk = maps_disk[image + '_' + ax_los]
if image == "Q":
im = P.imshow(
map_Q,
if image == 'Q' :
im = P.imshow(map_disk,
extent=im_extent,
origin="lower",
cmap="RdBu",
origin='lower',
cmap='RdBu',
norm=mpl.colors.LogNorm(),
vmin=0.01,
vmax=100.0,
)
vmin=0.01, vmax=100.)
else:
im = P.imshow(
map_disk,
@@ -348,47 +346,38 @@ def plot_maps(
cbar = P.colorbar(im)
if image == "coldens":
cbar.set_label(r"$log(N)$ (code)")
if image == 'coldens':
cbar.set_label(r'$log(\Sigma)$ (code)')
if "levels" in images:
map_level = maps_disk["levels_" + ax_los]
if 'levels' in images:
map_level = maps_disk['levels_' + ax_los]
# Computing linewidths
levels_ar = np.arange(np.min(map_level), np.max(map_level) + 1)
lw = np.ones(levels_ar.size) * 2
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]
)
lw[levels_ar < lvl_th] = 1.0
lw[levels_ar >= lvl_th] = lw[levels_ar >= lvl_th]**(lvl_th - levels_ar[levels_ar >= lvl_th])
lw[levels_ar < lvl_th] = 1.
cont = P.contour(
map_level,
cont = P.contour(map_level,
extent=im_extent,
origin="lower",
colors="k",
origin='lower',
colors='k',
linewidths=lw,
levels=levels_ar,
)
levels=levels_ar)
cont.levels = cont.levels + 1
P.clabel(
cont,
P.clabel(cont,
levels_ar[levels_ar < lvl_th + 2][1::2],
inline=1,
fontsize=8.0,
fmt="%1d",
)
elif image == "rho":
cbar.set_label(r"$log(n)$ (code)")
inline=1, fontsize=8., fmt='%1d')
elif image == 'rho':
cbar.set_label(r'$log(n)$ (code)')
if "speed" in images:
dmap_vh = maps_disk["v" + ax_h + "_" + ax_los]
dmap_vv = maps_disk["v" + ax_v + "_" + ax_los]
if 'speed' in images:
dmap_vh = maps_disk['v' + ax_h + '_' + ax_los]
dmap_vv = maps_disk['v' + ax_v + '_' + ax_los]
map_vh_red = dmap_vh[
::vel_red, ::vel_red
] # take only a subset of velocities
map_vv_red = dmap_vv[::vel_red, ::vel_red]
map_vh_red = dmap_vh[::vel_red,::vel_red] # take only a subset of velocities
map_vv_red = dmap_vv[::vel_red,::vel_red]
nh = map_vh_red.shape[0]
nv = map_vv_red.shape[1]
vec_h = (
@@ -532,9 +521,10 @@ def disk_prop(
g_az = (pos[:, 0] * g[:, 1] - pos[:, 1] * g[:, 0]) / norm_pos
# Select cells that are actually in the disk, ie within the scale height
G = 1.0
cs = np.sqrt(cells["P"] / cells["rho"]) # sound velocity
height = cs * np.sqrt(rc ** 3 / (G * mass_star))
G = 1.
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"] > 0.01 * np.mean(cells["rho"]) # condition on density
mask_vel = abs(v_rad / v_az) < 1.0 # condition on speed
@@ -559,8 +549,6 @@ def disk_prop(
total_mass_disk = np.sum(rho_disk * dvol_disk)
total_mass = np.sum(cells["rho"] * dx ** 3)
print("Mass disk", total_mass_disk)
print("Mass box", total_mass)
# Initialize binned quantities
cs_rad = np.zeros(nb_bin - 1)
@@ -726,38 +714,33 @@ def plot_disk_prop(
for plot in plots:
title = "t=" + str(time)[0:5] + " (code)"
P.grid()
P.xlabel("disk radius")
if plot == "rho":
P.xscale("log")
P.yscale("log")
P.plot(rad, prop_disk["rho"], color="k", linewidth=2)
P.ylabel(r"$n \, (code)$")
elif plot == "T":
P.xscale("log")
P.yscale("log")
P.plot(rad, prop_disk["temp"], color="k", linewidth=2)
P.ylabel(r"$T \, (K)$")
elif plot == "V":
P.xscale("log")
P.yscale("symlog", linthreshy=0.01)
P.plot(rad, prop_disk["v_rad"], color="k", linewidth=2)
P.plot(rad, prop_disk["v_kepl"], color="b", linewidth=2)
P.plot(rad, abs(prop_disk["v_az"]), color="r", linewidth=2)
P.plot(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"
)
P.ylabel(r"$V \, (km s^{-1})$")
elif plot == "coldens":
P.xscale("log")
P.yscale("log")
P.plot(rad, prop_disk["coldens"], color="k", linewidth=2)
P.ylabel(r"$N\, (cm^{-2})$")
elif plot == "alpha":
alpha_rey_mean, alpha_grav_mean = (
prop_disk["alpha_rey_mean"],
prop_disk["alpha_grav_mean"],
)
P.xlabel('disk radius')
if plot == 'rho':
P.xscale('log')
P.yscale('log')
P.plot(rad, prop_disk['rho'], color='k', linewidth=2)
P.ylabel(r'$n \, (code)$')
elif plot == 'T':
P.xscale('log')
P.yscale('log')
P.plot(rad,prop_disk['temp'],color='k',linewidth=2)
P.ylabel(r'$T \, (K)$')
elif plot == 'V':
P.xscale('log')
P.yscale('symlog',linthreshy=0.01)
P.plot(rad, prop_disk['v_rad'], color='k', linewidth=2)
P.plot(rad, prop_disk['v_kepl'], color='b', linewidth=2)
P.plot(rad, abs(prop_disk['v_az']), color='r', linewidth=2)
P.plot(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')
P.ylabel(r'$V \, (km s^{-1})$')
elif plot == 'coldens':
P.xscale('log')
P.yscale('log')
P.plot(rad,prop_disk['coldens'],color='k',linewidth=2)
P.ylabel(r'$\Sigma\, (cm^{-2})$')
elif plot == 'alpha':
alpha_rey_mean, alpha_grav_mean = prop_disk['alpha_rey_mean'], prop_disk['alpha_grav_mean']
P.xlim([1e-2, 0.25])
P.yscale("log")
P.ylim([1e-7, 1.0])
@@ -950,10 +933,10 @@ def disk_pdf(
nb_cells = np.sum(mask_flat)
P.grid()
P.yscale("log")
P.ylim([0.5 / nb_cells, 1.0])
P.xlabel(r"$\log(N / \bar{N})$")
P.ylabel(r"$\mathcal{P}_N$")
P.yscale('log')
P.ylim([0.5 / nb_cells, 1.])
P.xlabel(r'$\log(\Sigma / \bar{\Sigma})$')
P.ylabel(r'$\mathcal{P}_\Sigma$')
if put_title:
P.title(title)
values, edges, _ = P.hist(
@@ -999,14 +982,15 @@ def disk_pdf(
P.close()
# Derived quantities
drho = fluct_maps["rho"].flatten()
dcs = fluct_maps["cs"].flatten()
dv = fluct_maps["v"].flatten()
dvaz_kepl = abs(maps["vaz"] - v_kepl) / v_kepl
fluct_maps["vaz_kepl"] = dvaz_kepl
dmach = abs(maps["v"] - avg_maps["v"]) / maps["cs"]
fluct_maps["mach"] = dmach
dmach_mean = np.mean(dmach[mask_map].flatten())
drho = fluct_maps['rho'].flatten()
dcs = fluct_maps['cs'].flatten()
dv = fluct_maps['v'].flatten()
dvaz_kepl = abs(maps['vaz'] - v_kepl) / v_kepl
fluct_maps['vaz_kepl'] = dvaz_kepl
dmach = abs(maps['v'] - avg_maps['v']) / maps['cs']
dmach[dmach > 10.] = 10.
fluct_maps['mach'] = dmach
dmach_mean = np.mean(dmach[mask_map & (dmach < 5.)].flatten())
# Fluctuations plots
f = open(name_prop, "w")
@@ -1056,52 +1040,51 @@ def disk_pdf(
for cur_map in fluct_maps:
fluct_map = fluct_maps[cur_map]
if cur_map == "coldens":
im = P.imshow(
fluct_map,
if cur_map == 'coldens':
im = P.imshow(fluct_map,
norm=mpl.colors.LogNorm(),
extent=im_extent,
origin="lower",
cmap="viridis",
)
label = r"$log(N/\bar{N})$"
elif cur_map == "cs":
im = P.imshow(
np.log10(fluct_map),
origin='lower',
cmap='viridis')
label = r'$log(\Sigma/\bar{\Sigma})$'
elif cur_map == 'cs':
im = P.imshow(np.log10(fluct_map),
extent=im_extent,
origin="lower",
cmap="RdBu_r",
origin='lower',
cmap='RdBu_r',
vmin=-0.6,
vmax=0.6,
)
label = r"$log(c_s/\bar{c_s})$"
elif cur_map == "rho":
im = P.imshow(
np.log10(fluct_map),
vmax=0.6)
label = r'$log(c_s/\bar{c_s})$'
elif cur_map == 'rho':
im = P.imshow(np.log10(fluct_map),
extent=im_extent,
origin="lower",
cmap="RdBu_r",
vmin=-2.0,
vmax=2.0,
)
label = r"$log(\rho/\bar{\rho})$"
elif cur_map == "vaz_kepl":
im = P.imshow(
fluct_map,
origin='lower',
cmap='RdBu_r',
vmin=-2.,
vmax=2.)
label = r'$log(\rho/\bar{\rho})$'
elif cur_map == 'vaz_kepl':
im = P.imshow(fluct_map,
extent=im_extent,
origin="lower",
cmap="RdBu_r",
origin='lower',
cmap='RdBu_r',
norm=mpl.colors.LogNorm(),
vmax=1.0,
vmin=0.01,
)
label = r"$|v_\varphi - v_{kepl}|/v_{kepl}$"
elif cur_map == "vaz":
im = P.imshow(fluct_map, extent=im_extent, origin="lower", cmap="RdBu_r")
label = r"$v_\varphi / \bar{v_\varphi}$"
elif cur_map == "mach":
im = P.imshow(fluct_map, extent=im_extent, origin="lower", cmap="RdBu_r")
label = r"$|v - \bar{v}| / c_s$"
vmax=1.,
vmin=0.01)
label = r'$|v_\varphi - v_{kepl}|/v_{kepl}$'
elif cur_map == 'vaz':
im = P.imshow(fluct_map,
extent=im_extent,
origin='lower',
cmap='RdBu_r')
label = r'$v_\varphi / \bar{v_\varphi}$'
elif cur_map == 'mach':
im = P.imshow(fluct_map,
extent=im_extent,
origin='lower',
cmap='RdBu_r')
label = r'$|v - \bar{v}| / c_s$'
if put_title:
P.title(title)
@@ -1163,6 +1146,7 @@ def compare(
all_var = []
all_dmach = []
all_beta = []
for i, run in enumerate(runs):
path_run = path + "/" + run
@@ -1198,9 +1182,6 @@ def compare(
var_run.append(fit["var"])
dmach_run.append(prop_disk["dmach_mean"])
all_var = all_var + var_run
all_dmach = all_dmach + dmach_run
nb_outputs = nb_outputs + 1
print(run, num, nb_outputs)
@@ -1218,6 +1199,9 @@ def compare(
var_std[i] = np.std(var_run)
dmach[i] = np.mean(dmach_run)
dmach_std[i] = np.std(dmach_run)
all_var = all_var + var_run
all_dmach = all_dmach + dmach_run
all_beta = all_beta + [beta[i]]*len(var_run)
else:
for array in [alpha_rey, alpha_grav, Q]:
array[i] = np.nan
@@ -1229,10 +1213,10 @@ def compare(
if Q_in_name:
Q0[i] = float(run.split("_")[2][1:])
# Check if the output file exists, and exit if it is the case
name_save = path + "/alphaQ_" + nums_name + out_ext
# if (not force and len(glob.glob(name_save)) !=0):
# return
if pdf:
all_dmach = np.array(all_dmach)
all_var = np.array(all_var)
all_beta = np.array(all_beta)
# alpha = f(Qmin)
P.yscale("log")
@@ -1279,8 +1263,8 @@ def compare(
P.errorbar(beta, slope, yerr=slope_std, fmt="o")
P.legend()
P.ylabel(r"$d\log\mathcal{P}_N / d\logN$")
P.xlabel(r"$\beta$")
P.ylabel(r'$d\log\mathcal{P}_\Sigma / d\log\Sigma$')
P.xlabel(r'$\beta$')
(a, b, rho, _, stderr) = linregress(beta, slope)
P.plot(beta, a * beta + b, "--", linewidth=1.5)
@@ -1298,8 +1282,8 @@ def compare(
P.errorbar(beta, var, yerr=var_std, fmt="o")
P.legend()
P.ylabel(r"$Var(\log(N / \bar(N))$")
P.xlabel(r"$\beta$")
P.ylabel(r'$Var(\log(\Sigma / \bar{\Sigma})$')
P.xlabel(r'$\beta$')
(a, b, rho, _, stderr) = linregress(beta, var)
P.plot(beta, a * beta + b, "--", linewidth=1.5)
@@ -1313,17 +1297,33 @@ def compare(
P.close()
# var = f(log(dmach)
mask = all_dmach < 1
all_var = all_var[mask]
all_beta = all_beta[mask]
all_dmach = all_dmach[mask]
P.grid()
P.plot(all_dmach, all_var, "o")
cmap = mpl.cm.get_cmap('RdYlBu', np.max(beta) - np.min(beta) + 1)
P.scatter(all_dmach, np.exp(all_var), c=all_beta,
vmin=np.min(beta)-0.5,
vmax=np.max(beta)+0.5,
cmap=cmap)
cbar = P.colorbar()
cbar.set_ticks(beta)
cbar.set_label(r'$\beta$')
# P.errorbar(dmach, np.exp(var), xerr=dmach_std, yerr=np.exp(var) * var_std, fmt='+')
P.legend()
P.xlabel(r"$<(v - \bar{v}) / c_s>$")
P.ylabel(r"$Var(\log(N / \bar(N))$")
P.yscale("log")
P.xlabel(r'$<(v - \bar{v}) / c_s>$')
P.ylabel(r'$\exp(Var(\log(\Sigma / \bar{\Sigma}))$')
(a, b, rho, _, stderr) = linregress(all_dmach, np.exp(all_var))
#P.plot(all_dmach, a*all_dmach + b, '--', linewidth=1.5)
print("a=%e, b=%e, rho^2=%e"% (a,b,rho**2))
(a, b, rho, _, stderr) = linregress(dmach, np.exp(var))
# P.plot(all_dmach, a*all_dmach + b, '-.', linewidth=1.5)
print("a=%e, b=%e, rho^2=%e"% (a,b,rho**2))
# (a, b, rho, _, stderr) = linregress(var, log(dmach)
# P.plot(var, a*var + b, '--', linewidth=1.5)
# print("a=%e, b=%e, rho^2=%e"% (a,b,rho**2))
if interactive:
P.figure()
@@ -1467,11 +1467,12 @@ def evolution(path, nums, force=False, interactive=False, pdf=False):
print(time, slope)
P.plot(time, slope, "o-.")
P.legend()
P.ylabel(r"$d\log\mathcal{P}_{N} / d\logN$")
P.xlabel(r"time (code)")
P.ylabel(r'$d\log\mathcal{P}_{\Sigma} / d\log\Sigma$')
P.xlabel(r'time (code)')
if interactive:
P.figure()
else:
P.tight_layout(pad=1)
P.savefig(path + "/dcolslope_time" + out_ext)
P.close()
pipeline_disk.py
+1 -1
View File
@@ -246,7 +246,7 @@ for run in runs:
# If we are here, success !
success = True
run_succeded[run].append(i)
except (ValueError, IOError) as e:
except (ValueError, IOError, KeyError) as e:
print(e)
if args.watch and failures < args.allowed_failures:
failures = failures + 1