From 3afc117e2a11dc04cdb68680bc8911f94a55f82d Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Fri, 14 Jun 2019 14:14:02 +0200 Subject: [PATCH] Corrected errors --- disk_postprocess.py | 329 ++++++++++++++++++++++---------------------- pipeline_disk.py | 2 +- 2 files changed, 166 insertions(+), 165 deletions(-) diff --git a/disk_postprocess.py b/disk_postprocess.py index 07a801e..05ccaf6 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -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, - extent=im_extent, - origin="lower", - cmap="RdBu", - norm=mpl.colors.LogNorm(), - vmin=0.01, - vmax=100.0, - ) + if image == 'Q' : + im = P.imshow(map_disk, + extent=im_extent, + origin='lower', + cmap='RdBu', + norm=mpl.colors.LogNorm(), + 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] - # Computing linewidths - 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 + 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. - cont = P.contour( - map_level, - extent=im_extent, - origin="lower", - colors="k", - linewidths=lw, - levels=levels_ar, - ) - cont.levels = cont.levels + 1 + cont = P.contour(map_level, + extent=im_extent, + origin='lower', + colors='k', + linewidths=lw, + levels=levels_ar) + cont.levels = cont.levels + 1 - 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)") + P.clabel(cont, + levels_ar[levels_ar < lvl_th + 2][1::2], + 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,10 +521,11 @@ 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)) - mask_pos = np.abs(pos[:, 2]) < height # condition on position + 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, - 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), - extent=im_extent, - 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), - 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, - extent=im_extent, - 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$" + if cur_map == 'coldens': + im = P.imshow(fluct_map, + norm=mpl.colors.LogNorm(), + extent=im_extent, + 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', + 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), + extent=im_extent, + 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', + norm=mpl.colors.LogNorm(), + 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() + diff --git a/pipeline_disk.py b/pipeline_disk.py index b231186..0d8d457 100644 --- a/pipeline_disk.py +++ b/pipeline_disk.py @@ -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