From e654c5591dd0bed4089f23b84aa25d449907b168 Mon Sep 17 00:00:00 2001 From: Noe Brucy Date: Thu, 25 Apr 2019 09:31:21 +0200 Subject: [PATCH] Corrected error (lbox) in mass determination --- disk_postprocess.py | 81 +++++++++++++++++++++++++++++++++------------ 1 file changed, 59 insertions(+), 22 deletions(-) diff --git a/disk_postprocess.py b/disk_postprocess.py index 4b1e960..9a0592c 100644 --- a/disk_postprocess.py +++ b/disk_postprocess.py @@ -15,7 +15,7 @@ import glob as glob import pickle as pickle from pymses.sources.ramses import output -from pymses.analysis import Camera, raytracing, slicing +from pymses.analysis import Camera, raytracing, slicing, splatting from pymses.filters import CellsToPoints from pymses.analysis import ScalarOperator, FractionOperator, MaxLevelOperator @@ -39,6 +39,7 @@ def make_image_disk( cpuamr=False, pos_star=np.array([1.0, 1.0, 1.0]), interactive=False, + fft=False, ): """ Make several useful image of an output of a simulation @@ -77,6 +78,7 @@ def make_image_disk( put_title=put_title, pos_star=pos_star, interactive=interactive, + fft=fft, ) @@ -96,6 +98,7 @@ def make_image_aux( pos_star=np.array([1.0, 1.0, 1.0]), put_title=True, interactive=False, + fft=False, ): """ Make several useful image of an output of a simulation, auxillary function @@ -136,7 +139,12 @@ def make_image_aux( return rho_op = ScalarOperator(lambda dset: dset["rho"], ro.info["unit_density"]) - rt = raytracing.RayTracer(amr, ro.info, rho_op) + + rt = None + if fft: + rt = splatting.SplatterProcessor(amr, ro.info, rho_op) + else: + rt = raytracing.RayTracer(amr, ro.info, rho_op) axes_los = ["x", "y", "z"] # Line of sight axes axes_h = ["y", "x", "x"] # Horizontal axes @@ -336,7 +344,11 @@ def make_image_aux( ro.info["unit_velocity"], ) rt_omega = raytracing.RayTracer(amr, ro.info, omega_op) - rt_cs = raytracing.RayTracer(amr, ro.info, cs_op) + + if fft: + rt_cs = splatting.SplatterProcessor(amr, ro.info, cs_op, surf_qty=False) + else: + rt_cs = raytracing.RayTracer(amr, ro.info, cs_op) dmap_omega = rt_omega.process(cam) dmap_cs = rt_cs.process(cam) @@ -505,7 +517,7 @@ def disk_prop( amr = ro.amr_source(["rho", "vel", "Br", "Bl", "P"]) cell_source = CellsToPoints(amr) cells = cell_source.flatten() - dx = cells.get_sizes() + dx = cells.get_sizes() * lbox pos = cells.points * lbox # Get positions in the frame of the protostar pos = pos - pos_star @@ -540,8 +552,13 @@ def disk_prop( v_rad_disk = v_rad[mask] v_az_disk = v_az[mask] v_kepl = np.sqrt(mass_star * G / rc_disk) + height_disk = height[mask] 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) @@ -554,15 +571,12 @@ def disk_prop( v_rad_rad = np.zeros(nb_bin - 1) alpha_rey_rad = np.zeros(nb_bin - 1) Q_kepl_rad = np.zeros(nb_bin - 1) + height_rad = np.zeros(nb_bin - 1) for i in range(nb_bin - 1): mask_bin = (rc_disk > rad[i]) & (rc_disk < rad[i + 1]) - print( - "Bin #{} : {} cells between {} and {}".format( - i, np.sum(mask_bin), rad[i], rad[i + 1] - ) - ) + # print("Bin #{} : {} cells between {} and {}".format(i, 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] @@ -577,10 +591,8 @@ def disk_prop( # TODO verifier unites # Surface of a bin : S = dr * 2 * pi * r with # 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) ** 3 - / ((rad[i + 1] - rad[i]) * (rad[i + 1] + rad[i]) * np.pi) + coldens_rad[i] = np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) / ( + (rad[i + 1] - rad[i]) * (rad[i + 1] + rad[i]) * np.pi ) v_az_rad[i] = np.sum( @@ -590,6 +602,9 @@ def disk_prop( v_rad_rad[i] = np.sum( v_rad_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin] ) / np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) + height_rad[i] = np.sum( + height_disk[mask_bin] * rho_disk[mask_bin] * dvol_disk[mask_bin] + ) / np.sum(rho_disk[mask_bin] * dvol_disk[mask_bin]) alpha_rey_rad[i] = ( ( @@ -624,6 +639,7 @@ def disk_prop( prop_disk = { "time": time, "tot_mass": total_mass_disk, + "mass_box": total_mass, "rad": rad[0 : nb_bin - 1], "center": pos_star, "alpha_rey": alpha_rey_rad, @@ -636,6 +652,7 @@ def disk_prop( "temp": temp_rad, "cs": cs_rad, "Q_kepl": Q_kepl_rad, + "height": height_rad, } # store the results @@ -678,8 +695,10 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False): P.close() P.xscale("log") - P.plot(prop_disk["rad"], np.log10(prop_disk["rho"]), color="k", linewidth=2) - P.ylabel(r"$\log(n) \, (code)$") + P.yscale("log") + P.grid() + P.plot(prop_disk["rad"], prop_disk["rho"], color="k", linewidth=2) + P.ylabel(r"$n \, (code)$") P.xlabel("disk radius") P.title(title) if interactive: @@ -689,8 +708,10 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False): P.close() P.xscale("log") - P.plot(prop_disk["rad"], np.log10(prop_disk["temp"]), color="k", linewidth=2) - P.ylabel(r"$\log(T) \, (K)$") + P.yscale("log") + P.grid() + P.plot(prop_disk["rad"], prop_disk["temp"], color="k", linewidth=2) + P.ylabel(r"$T \, (K)$") P.xlabel("disk radius") P.title(title) if interactive: @@ -718,8 +739,10 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False): P.close() P.xscale("log") - P.plot(prop_disk["rad"], np.log10(prop_disk["coldens"]), color="k", linewidth=2) - P.ylabel(r"$\log(N) \, (cm^{-2})$") + P.yscale("log") + P.grid() + P.plot(prop_disk["rad"], prop_disk["coldens"], color="k", linewidth=2) + P.ylabel(r"$N\, (cm^{-2})$") P.xlabel("disk radius ") P.title(title) if interactive: @@ -763,6 +786,20 @@ def plot_disk_prop(path, num, force=False, tag="", interactive=False): P.savefig(path + "/Q_r_" + str(num).zfill(5) + out_ext) P.close() - if pdf: - P.savefig(path + "/Q_r_" + str(num).zfill(5) + ".pdf") - P.savefig(path + "/Q_r_" + str(num).zfill(5) + ".jpeg") + # height ration + P.grid() + P.plot( + prop_disk["rad"], + abs(prop_disk["height"] / prop_disk["rad"]), + color="b", + linewidth=2, + ) + P.ylabel(r"H ratio") + P.xlabel("disk radius ") + P.title(title + ", mass of box = {} (code)".format(prop_disk["mass_box"])) + + if interactive: + pass + else: + P.savefig(path + "/H_r_" + str(num).zfill(5) + out_ext) + P.close()