From aa7c541c25c1203af67033baffa89e5301bb462b Mon Sep 17 00:00:00 2001 From: Lukas Winkler Date: Tue, 23 Aug 2022 17:29:44 +0200 Subject: [PATCH] better comparisons of gas simulations --- auriga_comparison.py | 34 ++++++++++++++++--- slices.py | 78 +++++++++++++++++++------------------------- temperatures.py | 7 ++-- 3 files changed, 67 insertions(+), 52 deletions(-) diff --git a/auriga_comparison.py b/auriga_comparison.py index 6ee910f..c6ed13e 100644 --- a/auriga_comparison.py +++ b/auriga_comparison.py @@ -13,7 +13,7 @@ from matplotlib.axes import Axes from matplotlib.colors import LogNorm from matplotlib.figure import Figure -from cic import cic_from_radius +from cic import cic_from_radius, cic_range from halo_mass_profile import halo_mass_profile from nfw import fit_nfw from paths import auriga_dir, richings_dir @@ -168,9 +168,6 @@ for dir in sorted(root_dir.glob("*")): # linestyle="dotted" # ) - if has_baryons: - create_2d_slice(input_file, center, property="InternalEnergies") - centers[dir.name] = center if is_by_adrian: with reference_file.open("wb") as f: @@ -230,6 +227,35 @@ for dir in sorted(root_dir.glob("*")): ) ) i += 1 + + if has_baryons: + fig3, axs_baryon = plt.subplots(nrows=1, ncols=5, sharex="all", sharey="all", figsize=(10, 4)) + extent = [46, 52, 54, 60] # xrange[0], xrange[-1], yrange[0], yrange[-1] + for ii, property in enumerate(["cic", "Densities", "Entropies", "InternalEnergies", "Temperatures"]): + print(property) + if property == "cic": + grid, _ = cic_range(X + center[0], Y + center[1], 1000, *extent, periodic=False) + else: + grid = create_2d_slice(input_file, center, property=property, extent=extent) + print("minmax", grid.min(), grid.max()) + assert grid.min() != grid.max() + ax_baryon: Axes = axs_baryon[ii] + img = ax_baryon.imshow( + grid, + norm=LogNorm(), + interpolation="none", + origin="lower", + extent=extent, + ) + ax_baryon.set_title(property) + # ax_baryon.set_xlabel("X") + # ax_baryon.set_ylabel("Y") + ax_baryon.set_aspect("equal") + fig3.suptitle(input_file.parent.stem) + fig3.tight_layout() + fig3.savefig(Path("~/tmp/slice.png").expanduser(), dpi=300) + plt.show() + # plot_cic( # rho, extent, # title=str(dir.name) diff --git a/slices.py b/slices.py index 2cea9dd..1843a6a 100644 --- a/slices.py +++ b/slices.py @@ -1,5 +1,5 @@ from pathlib import Path -from typing import List +from typing import List, Tuple import h5py import matplotlib.pyplot as plt @@ -11,56 +11,46 @@ from temperatures import calculate_T from utils import create_figure +def filter_3d( + coords: np.ndarray, data: np.ndarray, + extent: List[float] +) -> Tuple[np.ndarray, np.ndarray]: + filter = ( + (extent[0] < coords[::, 0]) & + (coords[::, 0] < extent[1]) & + + (extent[2] < coords[::, 1]) & + (coords[::, 1] < extent[3]) + ) + print("before", coords.shape) + data = data[filter] + coords = coords[filter] + + print("after", coords.shape) + return coords, data + + def create_2d_slice( - input_file: Path, center: List[float], property: str, axis="Z", thickness=3, method="nearest" -): - axis_names = ["X", "Y", "Z"] - cut_axis = axis_names.index(axis) - limits = { - "X": (46, 52), - "Y": (54, 60), - "Z": (center[cut_axis] - 10, center[cut_axis] + 10) - } + input_file: Path, center: List[float], extent, + property="InternalEnergies", method="nearest" +) -> np.ndarray: + cut_axis = 2 # Z with h5py.File(input_file) as f: pt0 = f["PartType0"] coords = pt0["Coordinates"][:] - energies = pt0["InternalEnergies"][:] - entropies = pt0["Entropies"][:] + data = pt0[property if property != "Temperatures" else "InternalEnergies"][:] - print((center[cut_axis] - thickness < coords[::, cut_axis]).shape) - # in_slice = (center[cut_axis] - thickness < coords[::, cut_axis]) & ( - # coords[::, cut_axis] < center[cut_axis] + thickness) - # print("got slice") - # coords_in_slice = coords[in_slice] - # data_in_slice = data[in_slice] - filter = ( - (limits["X"][0] < coords[::, 0]) & - (coords[::, 0] < limits["X"][1]) & + coords, data = filter_3d(coords, data, extent) + if property == "Temperatures": + print("calculating temperatures") + data = np.array([calculate_T(u) for u in data]) - (limits["Y"][0] < coords[::, 1]) & - (coords[::, 1] < limits["Y"][1]) & - - (limits["Z"][0] < coords[::, 2]) & - (coords[::, 2] < limits["Z"][1]) - ) - print("before", coords.shape) - energies = energies[filter] - entropies = entropies[filter] - coords = coords[filter] - - print("after", coords.shape) - print("calculating temperatures") - temperatures = np.array([calculate_T(u) for u in energies]) - - other_axis = {"X": ("Y", "Z"), "Y": ("X", "Z"), "Z": ("X", "Y")} - x_axis_label, y_axis_label = other_axis[axis] - x_axis = axis_names.index(x_axis_label) - y_axis = axis_names.index(y_axis_label) - xrange = np.linspace(coords[::, x_axis].min(), coords[::, x_axis].max(), 1000) - yrange = np.linspace(coords[::, y_axis].min(), coords[::, y_axis].max(), 1000) + xrange = np.linspace(extent[0],extent[1], 1000) + yrange = np.linspace(extent[2],extent[3], 1000) gx, gy, gz = np.meshgrid(xrange, yrange, center[cut_axis]) print("interpolating") - grid = griddata(coords, temperatures, (gx, gy, gz), method=method)[::, ::, 0] + grid = griddata(coords, data, (gx, gy, gz), method=method)[::, ::, 0] + return grid print(grid.shape) # stats, x_edge, y_edge, _ = binned_statistic_2d( # coords_in_slice[::, x_axis], @@ -85,6 +75,4 @@ def create_2d_slice( ax.set_aspect("equal") fig.colorbar(img, label="Temperatures") fig.tight_layout() - fig.savefig(Path("~/tmp/slice.png").expanduser(), dpi=300) plt.show() - exit() diff --git a/temperatures.py b/temperatures.py index d2a0e5c..973dde6 100644 --- a/temperatures.py +++ b/temperatures.py @@ -16,9 +16,6 @@ const_boltzmann_k_cgs = 1.380649e-16 const_proton_mass = const_proton_mass_cgs / UnitMass_in_cgs const_boltzmann_k = const_boltzmann_k_cgs / UnitMass_in_cgs / UnitLength_in_cgs ** 2 * (UnitTime_in_cgs ** 2) -print(const_proton_mass) -print(const_boltzmann_k) -print() @njit @@ -33,7 +30,11 @@ def calculate_T(u): else: return T_transition + if __name__ == "__main__": + print(const_proton_mass) + print(const_boltzmann_k) + print() internal_energies = [6.3726251e+02, 7.7903375e+02, 1.7425287e+04, 6.4113910e+04, 3.8831848e+04, 1.1073163e+03, 7.7394878e+03, 7.5230023e+04, 9.1036992e+04, 2.4060946e+00]