From cfd008bc3e212a3a2406765327e108273415dd2b Mon Sep 17 00:00:00 2001 From: Lukas Winkler Date: Wed, 31 Aug 2022 15:26:27 +0200 Subject: [PATCH] major improvements to auriga/richings plots and barion data --- auriga_comparison.py | 218 +++++++++++++++++++++++++------------------ halo_mass_profile.py | 40 ++++---- paths.example.py | 2 + pbh_comparison.py | 11 +-- ramses.py | 11 ++- slices.py | 103 ++++++++++---------- utils.py | 3 + 7 files changed, 217 insertions(+), 171 deletions(-) diff --git a/auriga_comparison.py b/auriga_comparison.py index 3ab26a5..0132d79 100644 --- a/auriga_comparison.py +++ b/auriga_comparison.py @@ -4,6 +4,7 @@ from enum import Enum from pathlib import Path from pprint import pprint from subprocess import run +from sys import argv from typing import List, Tuple import h5py @@ -28,7 +29,7 @@ from nfw import fit_nfw from paths import auriga_dir, richings_dir from ramses import load_ramses_data, get_slice_argument, load_slice_data from readfiles import read_file, read_halo_file, ParticlesMeta -from slices import create_2d_slice +from slices import create_2d_slice, filter_3d from utils import read_swift_config, print_wall_time, figsize_from_page_fraction @@ -37,14 +38,22 @@ class Mode(Enum): auriga6 = 2 +class Plot(Enum): + auriga_plots = "auriga" + richings_bary = "richings_bary" + + mode = Mode.richings +try: + plottype = Plot(argv[1]) +except KeyError: + plottype = None + cache = HDFCache(Path("auriga_cache.hdf5")) -if mode == Mode.richings: - boxsize = 67.77 -else: - boxsize = None # not yet needed +if plottype == Plot.auriga_plots: + mode = Mode.auriga6 def dir_name_to_parameter(dir_name: str): @@ -68,7 +77,12 @@ def main(): ax1: Axes = fig1.gca() fig2: Figure = plt.figure(figsize=figsize_from_page_fraction()) ax2: Axes = fig2.gca() - fig4, axs_baryon = plt.subplots(nrows=2, ncols=4, sharex="all", sharey="all", figsize=(10, 4)) + axs_baryon: List[List[Axes]] + fig4, axs_baryon = plt.subplots( + nrows=2, ncols=4, + sharex="all", sharey="all", + figsize=figsize_from_page_fraction(columns=2, height_to_width=0.5) + ) fig5: Figure = plt.figure(figsize=figsize_from_page_fraction()) ax5: Axes = fig5.gca() fig6: Figure = plt.figure(figsize=figsize_from_page_fraction()) @@ -108,14 +122,18 @@ def main(): if not is_by_adrian: levelmin, levelmin_TF, levelmax = dir_name_to_parameter(dir.name) print(levelmin, levelmin_TF, levelmax) - if not has_baryons: - continue - if levelmax != 11: - continue - if not is_ramses: - continue + if plottype == Plot.auriga_plots: + if (levelmin, levelmin_TF, levelmax) == (7, 9, 9): + continue + elif plottype == Plot.richings_bary: + if not has_baryons: + continue + if levelmax != 11: + continue + # if not is_ramses: + # continue - input_file = dir / "output_0009.hdf5" + input_file = dir / "output_0007.hdf5" if mode == Mode.richings: input_file = dir / "output_0004.hdf5" if is_by_adrian or is_ramses: @@ -172,7 +190,8 @@ def main(): center = np.array([halo.X, halo.Y, halo.Z]) center = find_center(df, center) log_radial_bins, bin_masses, bin_densities, center = halo_mass_profile( - df, center, particles_meta, plot=False, num_bins=100, rmin=0.002, rmax=6.5 + df[["X", "Y", "Z"]].to_numpy(), center, particles_meta, plot=False, + num_bins=100, rmin=0.002, rmax=6.5 ) i_min_border = np.argmax( 0.01 < log_radial_bins @@ -195,9 +214,9 @@ def main(): with reference_file.open("wb") as f: pickle.dump([log_radial_bins, bin_masses, bin_densities], f) if is_by_adrian: - label = "reference" + label = "Reference" else: - label = f"{levelmin}, {levelmin_TF}, {levelmax}" + label = f"({levelmin}, {levelmin_TF}, {levelmax})" ax1.loglog(log_radial_bins[:-1], bin_masses, label=label, c=f"C{i}") ax2.loglog(log_radial_bins[:-1], bin_densities, label=label, c=f"C{i}") @@ -222,8 +241,7 @@ def main(): ax.axvline(4 * softening_length, color=f"C{i}", linestyle="dotted") # for ax in [ax1, ax2]: # ax.axvline(vr_halo.Rvir, color=f"C{i}", linestyle="dashed") - - X, Y, Z = df.X.to_numpy(), df.Y.to_numpy(), df.Z.to_numpy() + coords = df[["X", "Y", "Z"]].to_numpy() # shift: (-6, 0, -12) # if not is_by_adrian: @@ -232,28 +250,54 @@ def main(): # zshift = Zc - Zc_adrian # print("shift", xshift, yshift, zshift) - X -= center[0] - Y -= center[1] - Z -= center[2] + coords_centered = coords - center - rho, extent = cic_from_radius(X, Z, 4000, 0, 0, 5, periodic=False) + rho, extent = cic_from_radius(coords_centered[::, 0], coords_centered[::, 2], 500, 0, 0, 1.5, periodic=False) vmin = min(vmin, rho.min()) vmax = max(vmax, rho.max()) - - images.append( - Result( - rho=rho, - title=str(dir.name), - levels=(levelmin, levelmin_TF, levelmax) if levelmin else None, - ) + res = Result( + rho=rho, + title=f"levelmin={levelmin}, levelmin_TF={levelmin_TF}, levelmax={levelmax}" if not is_by_adrian else "Reference", + levels=(levelmin, levelmin_TF, levelmax) if not is_by_adrian else (100, 100, 100), ) + images.append(res) i += 1 if has_baryons: interpolation_method = "nearest" # "linear" bary_file = dir / "output_00009" if is_ramses else input_file - radius = 10 + if is_ramses: + s: RamsesSnap = pynbody.load(str(bary_file)) + gas_data: FamilySubSnap = s.gas + temperature_array: SimArray = gas_data["temp"] + p_array: SimArray = gas_data["p"].in_units("1e10 Msol Mpc^-3 km^2 s^-2") + rho_array: SimArray = gas_data["rho"].in_units("1e10 Msol Mpc^-3") + coord_array: SimArray = gas_data["pos"].in_units("Mpc") + mass_array = np.asarray(gas_data["mass"].in_units("1e10 Msol")) + bary_coords = np.asarray(coord_array) + bary_properties = { + "Temperatures": np.asarray(temperature_array.in_units("K")), + "Pressures": np.asarray(p_array), + "Densities": np.asarray(rho_array), + "Entropies": np.asarray(log10(p_array / rho_array ** (5 / 3))), + } + else: + with h5py.File(input_file) as f: + pt0 = f["PartType0"] + bary_coords = pt0["Coordinates"][:] + mass_array = pt0["Masses"][:] + bary_properties = { + "InternalEnergies": pt0["InternalEnergies"][:], + "Densities": pt0["Densities"][:], + "Pressures": pt0["Pressures"][:], + # "Entropies": log10(pt0["Densities"][:] / pt0["Densities"][:] ** (5 / 3)), + "Entropies": pt0["Entropies"][:] + } + bary_properties["Temperatures"] = bary_properties["InternalEnergies"] + + radius = 1.9 + resolution = 1000 # xrange[0], xrange[-1], yrange[0], yrange[-1] extent = [center[0] - radius, center[0] + radius, center[1] - radius, center[1] + radius] @@ -261,18 +305,22 @@ def main(): ramses_done = False for ii, property in enumerate(["cic", "Densities", "Entropies", "Temperatures"]): print("property:", property) - key = f"grid_{property}_{interpolation_method}_{radius}" + key = f"grid_{resolution}_{property}_{interpolation_method}_{radius}" cached_grid = cache.get(key, str(bary_file)) - if cached_grid is not None and not is_ramses: + if cached_grid is not None: grid = cached_grid else: + print("grid not yet cached, calculating now") if property == "cic": - rho, _ = cic_range(X + center[0], Y + center[1], 1000, *extent, periodic=False) - # TODO: filter in Z-axis + coords_in_box = filter_3d(coords, extent, zlimit=(center[2] - .1, center[2] + .1)) + rho, _ = cic_range(coords_in_box[::, 0], coords_in_box[::, 1], resolution, *extent, periodic=False) grid = 1.1 + rho.T else: if not is_ramses: - grid = create_2d_slice(bary_file, center, property=property, + grid = create_2d_slice(center, coords=bary_coords, + resolution=resolution, + property_name=property, + property_data=bary_properties[property], extent=extent, method=interpolation_method) else: frac_center = center / 100 @@ -280,24 +328,24 @@ def main(): print(frac_extent) print(frac_center) - args, imager_dir = get_slice_argument(frac_extent, frac_center, bary_file, - depth=.001) + args, imager_dir = get_slice_argument( + frac_extent, frac_center, + bary_file, depth=.001 + ) print(" ".join(args)) if not ramses_done: run(args, cwd=imager_dir) ramses_done = True property_map = { "Densities": "rhomap", - "Entropies": "Smap", # TODO: check - "InternalEnergies": None, + "Entropies": "Smap", "Temperatures": "Tmap" - } fname = imager_dir / f"snapshot_{property_map[property]}_zproj_zobs-0p00.bin" grid = load_slice_data(fname).T - cache.set(key, grid, str(input_file), compressed=True) - ax_baryon: Axes = axs_baryon[baryon_plot_counter, ii] + cache.set(key, grid, str(bary_file), compressed=True) + ax_baryon = axs_baryon[baryon_plot_counter][ii] img: AxesImage = ax_baryon.imshow( grid, norm=LogNorm(), @@ -305,56 +353,30 @@ def main(): origin="lower", extent=extent, ) - ax_baryon.set_title(property) + if baryon_plot_counter == 0: + ax_baryon.set_title(property) # ax_baryon.set_xlabel("X") # ax_baryon.set_ylabel("Y") ax_baryon.set_aspect("equal") # exit() baryon_plot_counter += 1 + continue - if not is_ramses: # TODO - continue - s: RamsesSnap = pynbody.load(str(bary_file)) - gas_data: FamilySubSnap = s.gas - temperature_array: SimArray = gas_data["temp"] - p_array: SimArray = gas_data["p"] - rho_array: SimArray = gas_data["rho"] - coord_array: SimArray = gas_data["pos"] - coordinates = np.asarray(coord_array.in_units("Mpc")) - properties = { - "temperature": np.asarray(temperature_array.in_units("K")), - "entropy": np.asarray(log10(p_array / rho_array ** (5 / 3))), - } + r, prof = property_profile(bary_coords, center, mass_array, bary_properties, num_bins=100, rmin=0.002, + rmax=6.5) + integrator_name = "Ramses" if is_ramses else "Swift" + label = f"{integrator_name} {levelmin}, {levelmin_TF}, {levelmax}" + ax5.set_title("Densities") + ax6.set_title("Pressures") + ax5.loglog(r[1:], prof["Densities"], label=label) + ax6.loglog(r[1:], prof["Pressures"], label=label) - r, prof = property_profile(coordinates, center, properties, num_bins=100, rmin=0.002, rmax=6.5) - ax5.loglog(r[1:], prof["temperature"]) - ax6.semilogx(r[1:], prof["entropy"]) - plt.show() - exit() - # # quick baryon profiles using pynbody - # s.gas["pos"] -= np.asarray(center) - # print("profile") - # p = Profile(s.gas, ndim=3) - # fig, ax = create_figure() - # ax5.plot(p['rbins'], p['density'], 'k') - # plt.show() - # exit() - - # plot_cic( - # rho, extent, - # title=str(dir.name) - # ) - ax1.legend() - ax2.legend() - fig1.tight_layout() - fig2.tight_layout() - - # fig3: Figure = plt.figure(figsize=(9, 9)) - # axes: List[Axes] = fig3.subplots(3, 3, sharex=True, sharey=True).flatten() fig3: Figure = plt.figure( - figsize=figsize_from_page_fraction(columns=2, height_to_width=1) + # just a bit more than 2/3 so that the two rows don't overlap + figsize=figsize_from_page_fraction(columns=2, height_to_width=33 / 48) ) - axes: List[Axes] = fig3.subplots(3, 3, sharex=True, sharey=True).flatten() + axes: List[Axes] = fig3.subplots(2, 3, sharex="all", sharey="all").flatten() + images.sort(key=lambda r: r.levels, reverse=True) for result, ax in zip(images, axes): data = 1.1 + result.rho @@ -365,21 +387,31 @@ def main(): norm=LogNorm(vmin=vmin_scaled, vmax=vmax_scaled), extent=extent, origin="lower", + cmap="Greys", interpolation="none" ) - ax.set_title(result.title) + ax.text( + 0.5, + 0.95, + result.title, + horizontalalignment="center", + verticalalignment="top", + transform=ax.transAxes, + ) - fig3.tight_layout() - fig3.subplots_adjust(right=0.825) - cbar_ax = fig3.add_axes([0.85, 0.05, 0.05, 0.9]) - fig3.colorbar(img, cax=cbar_ax) + for ax in [ax1, ax2, ax5, ax6]: + ax.legend() + for fig in [fig1, fig2, fig3, fig4, fig5, fig6]: + fig.tight_layout() + fig.subplots_adjust(wspace=0, hspace=0) + axs_baryon[0][0].set_ylabel("Swift") + axs_baryon[1][0].set_ylabel("Ramses") - fig1.savefig(Path(f"~/tmp/auriga1.pdf").expanduser()) - fig2.savefig(Path(f"~/tmp/auriga2.pdf").expanduser()) - fig3.savefig(Path("~/tmp/auriga3.pdf").expanduser()) + fig1.savefig(Path(f"~/tmp/{plottype.value}1.pdf").expanduser()) + fig2.savefig(Path(f"~/tmp/{plottype.value}2.pdf").expanduser()) + fig3.savefig(Path(f"~/tmp/{plottype.value}3.pdf").expanduser()) - fig4.tight_layout() - fig4.savefig(Path("~/tmp/slice.png").expanduser(), dpi=300) + fig4.savefig(Path(f"~/tmp/{plottype.value}4.pdf").expanduser()) pprint(centers) plt.show() diff --git a/halo_mass_profile.py b/halo_mass_profile.py index 9c1e8df..a7a77f2 100644 --- a/halo_mass_profile.py +++ b/halo_mass_profile.py @@ -4,12 +4,12 @@ from typing import Dict import matplotlib.pyplot as plt import numpy as np -import pandas as pd from matplotlib.axes import Axes from matplotlib.figure import Figure -from scipy.spatial import KDTree from readfiles import ParticlesMeta, read_file, read_halo_file +from temperatures import calculate_T +from utils import print_progress def V(r): @@ -17,7 +17,7 @@ def V(r): def halo_mass_profile( - particles: pd.DataFrame, + positions: np.ndarray, center: np.ndarray, particles_meta: ParticlesMeta, rmin: float, @@ -25,9 +25,7 @@ def halo_mass_profile( plot=False, num_bins=30, ): - positions = particles[["X", "Y", "Z"]].to_numpy() distances = np.linalg.norm(positions - center, axis=1) - group_radius = distances.max() log_radial_bins = np.geomspace(rmin, rmax, num_bins) @@ -65,27 +63,31 @@ def halo_mass_profile( return log_radial_bins, bin_masses, bin_densities, center -def property_profile(positions: np.ndarray, center: np.ndarray, properties: Dict[str, np.ndarray], +def property_profile(positions: np.ndarray, center: np.ndarray, masses: np.ndarray, properties: Dict[str, np.ndarray], rmin: float, rmax: float, num_bins: int): - print("building KDTree") - tree = KDTree(positions) - print("done") - + distances = np.linalg.norm(positions - center, axis=1) log_radial_bins = np.geomspace(rmin, rmax, num_bins) - particles_inner_ring = set(tree.query_ball_point(center, rmin)) means = {} for key in properties.keys(): means[key] = [] - for r in log_radial_bins[1:]: - print(r) - particles_inside = set(tree.query_ball_point(center, r)) - particles_in_ring = particles_inside - particles_inner_ring + for k in range(num_bins - 1): + bin_start = log_radial_bins[k] + bin_end = log_radial_bins[k + 1] + print_progress(k, num_bins - 2, bin_end) + in_bin = np.where((bin_start < distances) & (distances < bin_end))[0] + masses_in_ring = masses[in_bin] for property, property_values in properties.items(): - prop_in_ring = property_values[list(particles_in_ring)] - mean_in_ring = np.mean(prop_in_ring) + if property == "InternalEnergies": + continue + prop_in_ring = property_values[in_bin] + if property == "Temperatures": + prop_in_ring = np.array([calculate_T(u) for u in prop_in_ring]) + + # mean_in_ring_unweighted = np.mean(prop_in_ring) + mean_in_ring = (prop_in_ring * masses_in_ring).sum() / masses_in_ring.sum() + # print(mean_in_ring_unweighted, mean_in_ring) means[property].append(mean_in_ring) - particles_inner_ring = particles_inside return log_radial_bins, means @@ -105,4 +107,4 @@ if __name__ == "__main__": halo = df_halos.loc[halo_id] - halo_mass_profile(particles_in_halo, halo, particles_meta, plot=True) + halo_mass_profile(particles_in_halo[["X", "Y", "Z"]].to_numpy(), halo, particles_meta, plot=True) diff --git a/paths.example.py b/paths.example.py index 73434ec..0e23dd1 100644 --- a/paths.example.py +++ b/paths.example.py @@ -6,3 +6,5 @@ auriga_dir = Path("/path/to/auriga6/") richings_dir = Path("path/to/richings21_ics/").expanduser() spectra_dir = Path("/path/to/git/spectra/build/") has_1024_simulations = False +ramses_imager = Path("~/cosmoca/RamsesImager/ramses_imager").expanduser() +pbh_dir = Path("~/cosmos_data/PBH/").expanduser() diff --git a/pbh_comparison.py b/pbh_comparison.py index 33cd351..1f54835 100644 --- a/pbh_comparison.py +++ b/pbh_comparison.py @@ -1,18 +1,16 @@ import numpy as np -import pandas as pd from matplotlib import pyplot as plt from matplotlib.axes import Axes from matplotlib.figure import Figure from cic import cic_from_radius, plot_cic -from find_center import find_center from halo_mass_profile import halo_mass_profile from paths import pbh_dir from readfiles import read_g4_file, ParticlesMeta from utils import figsize_from_page_fraction -def cic_comparison(pbh_high_coord, ref_high_coord,center): +def cic_comparison(pbh_high_coord, ref_high_coord, center): rhos = [] i = 0 for coord in [ref_high_coord, pbh_high_coord]: @@ -38,7 +36,7 @@ def main(): pbh_dir / "10000sigma" / "snapshot_005.hdf5", zoom_type="pbh") center = [30, 32, 30] - cic_comparison(ref_data[0], pbh_data[0],center) + cic_comparison(ref_data[0], pbh_data[0], center) fig1: Figure = plt.figure(figsize=figsize_from_page_fraction()) ax1: Axes = fig1.gca() fig2: Figure = plt.figure(figsize=figsize_from_page_fraction()) @@ -46,7 +44,6 @@ def main(): centered = False for data in [ref_data, pbh_data]: highres_coords, lowres_coords, highres_mass, lowres_mass = data - df = pd.DataFrame(highres_coords, columns=["X", "Y", "Z"]) particles_meta = ParticlesMeta(particle_mass=highres_mass) # center = np.median(highres_coords, axis=0) print(center) @@ -54,13 +51,13 @@ def main(): # center = find_center(df, center, initial_radius=0.01) centered = True log_radial_bins, bin_masses, bin_densities, center = halo_mass_profile( - df, center, particles_meta, plot=False, num_bins=100, vmin=0.002, vmax=5 + highres_coords, center, particles_meta, plot=False, num_bins=100, rmin=0.002, rmax=5 ) ax1.loglog(log_radial_bins[:-1], bin_masses) ax2.loglog(log_radial_bins[:-1], bin_densities) plt.show() - cic_comparison(ref_data[0], pbh_data[0],center) + cic_comparison(ref_data[0], pbh_data[0], center) if __name__ == '__main__': diff --git a/ramses.py b/ramses.py index 31a1863..2a44ebc 100644 --- a/ramses.py +++ b/ramses.py @@ -1,12 +1,15 @@ from pathlib import Path from typing import List +import matplotlib.pyplot as plt import numpy as np import pynbody +from pynbody.analysis.profile import Profile from pynbody.array import SimArray from pynbody.snapshot.ramses import RamsesSnap from readfiles import ParticlesMeta +from utils import create_figure def load_ramses_data(ramses_dir: Path): @@ -15,7 +18,12 @@ def load_ramses_data(ramses_dir: Path): coord_array: SimArray = s.dm["pos"] a = s.properties["a"] print("RAMSES a", a) - + # p = Profile(s.gas, ndim=3) + # s.gas["pos"]-= + # fig,ax=create_figure() + # ax.plot(p['rbins'], p['density'], 'k') + # plt.show() + # exit() masses = np.asarray(mass_array.in_units("1e10 Msol")) high_res_mass = np.amin(np.unique(masses)) # get lowest mass of particles is_high_res_particle = masses == high_res_mass @@ -25,6 +33,7 @@ def load_ramses_data(ramses_dir: Path): particles_meta = ParticlesMeta(particle_mass=high_res_mass) center = np.median(hr_coordinates, axis=0) + return hr_coordinates, particles_meta, center diff --git a/slices.py b/slices.py index 1843a6a..eefc0fb 100644 --- a/slices.py +++ b/slices.py @@ -1,7 +1,5 @@ -from pathlib import Path from typing import List, Tuple -import h5py import matplotlib.pyplot as plt import numpy as np from matplotlib.colors import LogNorm @@ -12,9 +10,9 @@ from utils import create_figure def filter_3d( - coords: np.ndarray, data: np.ndarray, - extent: List[float] -) -> Tuple[np.ndarray, np.ndarray]: + coords: np.ndarray, extent: List[float], data: np.ndarray = None, zlimit=None + +) -> Tuple[np.ndarray, np.ndarray] | np.ndarray: filter = ( (extent[0] < coords[::, 0]) & (coords[::, 0] < extent[1]) & @@ -22,57 +20,60 @@ def filter_3d( (extent[2] < coords[::, 1]) & (coords[::, 1] < extent[3]) ) + if zlimit: + filter = filter & ( + (zlimit[0] < coords[::, 2]) & + (coords[::, 2] < zlimit[1]) + + ) print("before", coords.shape) - data = data[filter] + if data is not None: + data = data[filter] coords = coords[filter] print("after", coords.shape) - return coords, data + if data is not None: + return coords, data + return coords -def create_2d_slice( - input_file: Path, center: List[float], extent, - property="InternalEnergies", method="nearest" -) -> np.ndarray: +def create_2d_slice(center: List[float], extent, coords: np.ndarray, property_name: str, property_data: np.ndarray, + resolution: int, + method="nearest") -> np.ndarray: cut_axis = 2 # Z - with h5py.File(input_file) as f: - pt0 = f["PartType0"] - coords = pt0["Coordinates"][:] - data = pt0[property if property != "Temperatures" else "InternalEnergies"][:] - coords, data = filter_3d(coords, data, extent) - if property == "Temperatures": - print("calculating temperatures") - data = np.array([calculate_T(u) for u in data]) - - 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, 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], - # coords_in_slice[::, y_axis], - # data_in_slice, - # bins=500, - # statistic="mean" - # ) - fig, ax = create_figure() - # stats = np.nan_to_num(stats) - print("plotting") - img = ax.imshow( - grid, - norm=LogNorm(), - interpolation="nearest", - origin="lower", - extent=[xrange[0], xrange[-1], yrange[0], yrange[-1]], - ) - ax.set_title(input_file.parent.stem) - ax.set_xlabel(x_axis_label) - ax.set_ylabel(y_axis_label) - ax.set_aspect("equal") - fig.colorbar(img, label="Temperatures") - fig.tight_layout() - plt.show() + coords, property_data = filter_3d(coords, extent, property_data) + if property_name == "Temperatures": + print("calculating temperatures") + property_data = np.array([calculate_T(u) for u in property_data]) + xrange = np.linspace(extent[0], extent[1], resolution) + yrange = np.linspace(extent[2], extent[3], resolution) + gx, gy, gz = np.meshgrid(xrange, yrange, center[cut_axis]) + print("interpolating") + grid = griddata(coords, property_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], + # coords_in_slice[::, y_axis], + # data_in_slice, + # bins=500, + # statistic="mean" + # ) + fig, ax = create_figure() + # stats = np.nan_to_num(stats) + print("plotting") + img = ax.imshow( + grid, + norm=LogNorm(), + interpolation="nearest", + origin="lower", + extent=[xrange[0], xrange[-1], yrange[0], yrange[-1]], + ) + ax.set_title(input_file.parent.stem) + ax.set_xlabel(x_axis_label) + ax.set_ylabel(y_axis_label) + ax.set_aspect("equal") + fig.colorbar(img, label="Temperatures") + fig.tight_layout() + plt.show() diff --git a/utils.py b/utils.py index 2dd348f..6df8aaf 100644 --- a/utils.py +++ b/utils.py @@ -39,6 +39,9 @@ def read_swift_config(dir: Path): def print_wall_time(dir: Path): + """ + Attention: This idea is flawed as it only shows the wall time of the last time the simulation was restarted + """ with (dir / "swift.log").open() as f: last_line = f.readlines()[-1] print(last_line)