Commit ba958c91 authored by Josh Borrow's avatar Josh Borrow Committed by Matthieu Schaller
Browse files

Adds ANARCHY-PU Scheme

parent 72de43e0
......@@ -1239,7 +1239,7 @@ esac
# Hydro scheme.
AC_ARG_WITH([hydro],
[AS_HELP_STRING([--with-hydro=<scheme>],
[Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, debug default: gadget2@:>@]
[Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, anarchy-pu debug default: gadget2@:>@]
)],
[with_hydro="$withval"],
[with_hydro="gadget2"]
......@@ -1284,6 +1284,9 @@ case "$with_hydro" in
planetary)
AC_DEFINE([PLANETARY_SPH], [1], [Planetary SPH])
;;
anarchy-pu)
AC_DEFINE([ANARCHY_PU_SPH], [1], [ANARCHY (PU) SPH])
;;
*)
......
......@@ -236,7 +236,7 @@ implicit problem. A valid section of the YAML file looks like:
H_reion_z: 11.5 # Redhift of Hydrogen re-ionization
He_reion_z_centre: 3.5 # Centre of the Gaussian used for Helium re-ionization
He_reion_z_sigma: 0.5 # Width of the Gaussian used for Helium re-ionization
He_reion_eV_p_H: 2.0 # Energy injected in eV per Hydrogen atom for Helium re-ionization.
He_reion_ev_p_H: 2.0 # Energy injected in eV per Hydrogen atom for Helium re-ionization.
And the optional parameters are:
......
python3 makeImage.py santabarbara_0153.hdf5 0 twilight white
python3 plotSolution.py 153 halo
python3 plotTempEvolution.py
python3 rhoTHaloComparison.py
"""
Plots the smoothing length (compared to the softening).
"""
import numpy as np
import matplotlib.pyplot as plt
import h5py
from collections import namedtuple
SnapshotData = namedtuple(
"SnapshotData",
[
"smoothing_lengths",
"particle_ids",
"softening",
"internal_length",
"snapshot_length",
],
)
HaloCatalogueData = namedtuple(
"HaloCatalogueData", ["largest_halo", "particle_ids_in_largest_halo"]
)
def load_data(filename: str) -> SnapshotData:
"""
Loads the data that we need, i.e. the smoothing lengths and the
softening length, from the snapshot.
"""
with h5py.File(filename, "r") as handle:
data = SnapshotData(
smoothing_lengths=handle["PartType0/SmoothingLength"][...],
particle_ids=handle["PartType0/ParticleIDs"][...],
softening=handle["GravityScheme"].attrs[
"Comoving softening length [internal units]"
][0],
internal_length=handle["InternalCodeUnits"].attrs[
"Unit length in cgs (U_L)"
][0],
snapshot_length=handle["Units"].attrs["Unit length in cgs (U_L)"][0],
)
return data
def load_halo_data(halo_filename: str) -> HaloCatalogueData:
"""
Loads the halo data and finds the particle IDs that belong to
the largest halo. The halo filename should be given without
any extension as we need a couple of files to complete this.
"""
catalogue_filename = f"{halo_filename}.properties"
groups_filename = f"{halo_filename}.catalog_groups"
particles_filename = f"{halo_filename}.catalog_particles"
with h5py.File(catalogue_filename, "r") as handle:
largest_halo = np.where(
handle["Mass_200crit"][...] == handle["Mass_200crit"][...].max()
)[0][0]
with h5py.File(groups_filename, "r") as handle:
offset_begin = handle["Offset"][largest_halo]
offset_end = handle["Offset"][largest_halo + 1]
with h5py.File(particles_filename, "r") as handle:
particle_ids = handle["Particle_IDs"][offset_begin:offset_end]
return HaloCatalogueData(
largest_halo=largest_halo, particle_ids_in_largest_halo=particle_ids
)
def make_plot(
snapshot_filename: str,
halo_filename: str,
output_filename="smoothing_length_variation.png",
) -> None:
"""
Makes the plot and saves it in output_filename.
The halo filename should be provided without extension.
"""
data = load_data(filename)
halo_data = load_halo_data(halo_filename)
smoothing_lengths_in_halo = data.smoothing_lengths[
np.in1d(data.particle_ids, halo_data.particle_ids_in_largest_halo)
]
softening = data.softening * (data.snapshot_length / data.internal_length)
fig, ax = plt.subplots(1)
ax.semilogy()
ax.hist(data.smoothing_lengths, bins="auto", label="All particles")
ax.hist(
smoothing_lengths_in_halo,
bins="auto",
label=f"Particles in largest halo (ID={halo_data.largest_halo})",
)
ax.axvline(x=softening, label="Softening", ls="--", color="grey")
ax.legend()
ax.set_xlabel("Smoothing length")
ax.set_ylabel("Number of particles")
ax.set_xlim(0, ax.get_xlim()[1])
fig.tight_layout()
fig.savefig(output_filename, dpi=300)
return
if __name__ == "__main__":
import argparse as ap
PARSER = ap.ArgumentParser(
description="""
Makes a plot of the smoothing lengths in the box, compared
to the gravitational softening. Also splits out the particles
that are contained in the largest halo according to the
velociraptor outputs.
"""
)
PARSER.add_argument(
"-s",
"--snapshot",
help="""
Filename and path for the snapshot (without the .hdf5),
Default: ./santabarbara_0153
""",
required=False,
default="./santabarbara_0153",
)
PARSER.add_argument(
"-v",
"--velociraptor",
help="""
The filename and path of the velociraptor files, excluding
the descriptors (i.e. without .catalog_particles).
Default: ./halo/santabarbara
""",
required=False,
default="./halo/santabarbara",
)
ARGS = vars(PARSER.parse_args())
filename = f"{ARGS['snapshot']}.hdf5"
make_plot(filename, ARGS["velociraptor"])
......@@ -106,13 +106,18 @@ def get_halo_data(catalogue_filename: str) -> HaloData:
that is given by VELOCIraptor.
"""
with h5py.File(catalogue_filename, "r") as file:
x = file["Xc"][0]
y = file["Yc"][0]
z = file["Zc"][0]
Mvir = file["Mass_200crit"][0]
Rvir = file["R_200crit"][0]
c = file["cNFW"][0]
largest_halo = np.where(
file["Mass_200crit"][...] == file["Mass_200crit"][...].max()
)
x = float(np.take(file["Xc"], largest_halo))
y = float(np.take(file["Yc"], largest_halo))
z = float(np.take(file["Zc"], largest_halo))
Mvir = float(np.take(file["Mass_200crit"], largest_halo))
Rvir = float(np.take(file["R_200crit"], largest_halo))
c = float(np.take(file["cNFW"], largest_halo))
return HaloData(c=c, Rvir=Rvir, Mvir=Mvir, center=np.array([x, y, z]))
......
################################################################################
# This file is part of SWIFT.
# Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
# Computes the temperature evolution of the gas in a cosmological box
# Physical constants needed for internal energy to temperature conversion
k_in_J_K = 1.38064852e-23
mH_in_kg = 1.6737236e-27
# Number of snapshots generated
n_snapshots = 153
snapname = "santabarbara"
import matplotlib
matplotlib.use("Agg")
from pylab import *
import h5py
import os.path
# Plot parameters
params = {'axes.labelsize': 10,
'axes.titlesize': 10,
'font.size': 9,
'legend.fontsize': 9,
'xtick.labelsize': 10,
'ytick.labelsize': 10,
'text.usetex': True,
'figure.figsize' : (3.15,3.15),
'figure.subplot.left' : 0.14,
'figure.subplot.right' : 0.99,
'figure.subplot.bottom' : 0.12,
'figure.subplot.top' : 0.99,
'figure.subplot.wspace' : 0.15,
'figure.subplot.hspace' : 0.12,
'lines.markersize' : 6,
'lines.linewidth' : 2.,
'text.latex.unicode': True
}
rcParams.update(params)
rc('font',**{'family':'sans-serif','sans-serif':['Times']})
# Read the simulation data
sim = h5py.File("%s_0000.hdf5" % snapname, "r")
boxSize = sim["/Header"].attrs["BoxSize"][0]
time = sim["/Header"].attrs["Time"][0]
scheme = sim["/HydroScheme"].attrs["Scheme"][0]
kernel = sim["/HydroScheme"].attrs["Kernel function"][0]
neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
H_transition_temp = sim["/HydroScheme"].attrs["Hydrogen ionization transition temperature"][0]
T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
git = sim["Code"].attrs["Git Revision"]
# Cosmological parameters
H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
unit_length_in_si = 0.01 * unit_length_in_cgs
unit_mass_in_si = 0.001 * unit_mass_in_cgs
unit_time_in_si = unit_time_in_cgs
# Primoridal ean molecular weight as a function of temperature
def mu(T, H_frac=H_mass_fraction, T_trans=H_transition_temp):
if T > T_trans:
return 4. / (8. - 5. * (1. - H_frac))
else:
return 4. / (1. + 3. * H_frac)
# Temperature of some primoridal gas with a given internal energy
def T(u, H_frac=H_mass_fraction, T_trans=H_transition_temp):
T_over_mu = (gas_gamma - 1.) * u * mH_in_kg / k_in_J_K
ret = np.ones(np.size(u)) * T_trans
# Enough energy to be ionized?
mask_ionized = (T_over_mu > (T_trans+1) / mu(T_trans+1, H_frac, T_trans))
if np.sum(mask_ionized) > 0:
ret[mask_ionized] = T_over_mu[mask_ionized] * mu(T_trans*10, H_frac, T_trans)
# Neutral gas?
mask_neutral = (T_over_mu < (T_trans-1) / mu((T_trans-1), H_frac, T_trans))
if np.sum(mask_neutral) > 0:
ret[mask_neutral] = T_over_mu[mask_neutral] * mu(0, H_frac, T_trans)
return ret
z = np.zeros(n_snapshots)
a = np.zeros(n_snapshots)
T_mean = np.zeros(n_snapshots)
T_std = np.zeros(n_snapshots)
T_log_mean = np.zeros(n_snapshots)
T_log_std = np.zeros(n_snapshots)
T_median = np.zeros(n_snapshots)
T_min = np.zeros(n_snapshots)
T_max = np.zeros(n_snapshots)
# Loop over all the snapshots
for i in range(n_snapshots):
sim = h5py.File("%s_%04d.hdf5"% (snapname, i), "r")
z[i] = sim["/Cosmology"].attrs["Redshift"][0]
a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
u = sim["/PartType0/InternalEnergy"][:]
# Compute the temperature
u *= (unit_length_in_si**2 / unit_time_in_si**2)
u /= a[i]**(3 * (gas_gamma - 1.))
Temp = T(u)
# Gather statistics
T_median[i] = np.median(Temp)
T_mean[i] = Temp.mean()
T_std[i] = Temp.std()
T_log_mean[i] = np.log10(Temp).mean()
T_log_std[i] = np.log10(Temp).std()
T_min[i] = Temp.min()
T_max[i] = Temp.max()
# CMB evolution
a_evol = np.logspace(-3, 0, 60)
T_cmb = (1. / a_evol)**2 * 2.72
# Plot the interesting quantities
figure()
subplot(111, xscale="log", yscale="log")
fill_between(a, T_mean-T_std, T_mean+T_std, color='C0', alpha=0.1)
plot(a, T_max, ls='-.', color='C0', lw=1., label="${\\rm max}~T$")
plot(a, T_min, ls=':', color='C0', lw=1., label="${\\rm min}~T$")
plot(a, T_mean, color='C0', label="${\\rm mean}~T$", lw=1.5)
fill_between(a, 10**(T_log_mean-T_log_std), 10**(T_log_mean+T_log_std), color='C1', alpha=0.1)
plot(a, 10**T_log_mean, color='C1', label="${\\rm mean}~{\\rm log} T$", lw=1.5)
plot(a, T_median, color='C2', label="${\\rm median}~T$", lw=1.5)
legend(loc="upper left", frameon=False, handlelength=1.5)
# Expected lines
plot([1e-10, 1e10], [H_transition_temp, H_transition_temp], 'k--', lw=0.5, alpha=0.7)
text(2.5e-2, H_transition_temp*1.07, "$T_{\\rm HII\\rightarrow HI}$", va="bottom", alpha=0.7, fontsize=8)
plot([1e-10, 1e10], [T_minimal, T_minimal], 'k--', lw=0.5, alpha=0.7)
text(1e-2, T_minimal*0.8, "$T_{\\rm min}$", va="top", alpha=0.7, fontsize=8)
plot(a_evol, T_cmb, 'k--', lw=0.5, alpha=0.7)
text(a_evol[20], T_cmb[20]*0.55, "$(1+z)^2\\times T_{\\rm CMB,0}$", rotation=-34, alpha=0.7, fontsize=8, va="top", bbox=dict(facecolor='w', edgecolor='none', pad=1.0, alpha=0.9))
redshift_ticks = np.array([0., 1., 2., 5., 10., 20., 50., 100.])
redshift_labels = ["$0$", "$1$", "$2$", "$5$", "$10$", "$20$", "$50$", "$100$"]
a_ticks = 1. / (redshift_ticks + 1.)
xticks(a_ticks, redshift_labels)
minorticks_off()
xlabel("${\\rm Redshift}~z$", labelpad=0)
ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
xlim(9e-3, 1.1)
ylim(20, 2.5e7)
savefig("Temperature_evolution.png", dpi=200)
"""
This script finds the temperatures inside all of the halos and
compares it against the virial temperature. This uses velociraptor
and the SWIFT snapshot.
Folkert Nobels (2018) nobels@strw.leidenuniv.nl
Josh Borrow (2019) joshua.borrow@durham.ac.uk
"""
import numpy as np
import h5py
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
mH = 1.6733e-24 # g
kB = 1.38e-16 # erg/K
def virial_temp(mu, M, h=0.703, a=1.0):
"""
Calculates the virial temperature according to
https://arxiv.org/pdf/1105.5701.pdf
Equation 1.
"""
return 4e4 * (mu / 1.2) * (M * h / 1e8) ** (2 / 3) / (10 * a)
def calculate_group_sizes_array(offsets: np.array, total_size: int) -> np.array:
"""
Calculates the group sizes array from the offsets and total size, i.e. it
calculates the diff between all of the offsets.
"""
# Does not include the LAST one
group_sizes = [x - y for x, y in zip(offsets[1:], offsets[:-1])]
group_sizes += [total_size - offsets[-1]]
group_sizes = np.array(group_sizes, dtype=type(offsets[0]))
return group_sizes
def create_group_array(group_sizes: np.array) -> np.array:
"""
Creates an array that looks like:
[GroupID0, GroupID0, ..., GroupIDN, GroupIDN]
i.e. for each group create the correct number of group ids.
This is used to be sorted alongside the particle IDs to track
the placement of group IDs.
"""
slices = []
running_total_of_particles = type(offsets[0])(0)
for group in group_sizes:
slices.append([running_total_of_particles, group + running_total_of_particles])
running_total_of_particles += group
groups = np.empty(group_sizes.sum(), dtype=int)
for group_id, group in enumerate(slices):
groups[group[0] : group[1]] = group_id
return groups
if __name__ == "__main__":
import argparse as ap
PARSER = ap.ArgumentParser(
description="""
Makes many plots comparing the virial temperature and the
temperature of halos. Requires the velociraptor files and
the SWIFT snapshot.
"""
)
PARSER.add_argument(
"-s",
"--snapshot",
help="""
Filename and path for the snapshot (without the .hdf5),
Default: ./santabarbara_0153
""",
required=False,
default="./santabarbara_0153",
)
PARSER.add_argument(
"-v",
"--velociraptor",
help="""
The filename and path of the velociraptor files, excluding
the descriptors (i.e. without .catalog_particles).
Default: ./halo/santabarbara
""",
required=False,
default="./halo/santabarbara",
)
ARGS = vars(PARSER.parse_args())
# Grab some metadata before we begin.
with h5py.File("%s.hdf5" % ARGS["snapshot"], "r") as handle:
# Cosmology
a = handle["Cosmology"].attrs["Scale-factor"][0]
h = handle["Cosmology"].attrs["h"][0]
# Gas
hydro = handle["HydroScheme"].attrs
X = hydro["Hydrogen mass fraction"][0]
gamma = hydro["Adiabatic index"][0]
mu = 1 / (X + (1 - X) / 4)
# First we must construct a group array so we know which particles belong
# to which group.
with h5py.File("%s.catalog_groups" % ARGS["velociraptor"], "r") as handle:
offsets = handle["Offset"][...]
# Then, extract the particles that belong to the halos. For that, we need
# the particle IDs:
with h5py.File("%s.catalog_particles" % ARGS["velociraptor"], "r") as handle:
ids_in_halos = handle["Particle_IDs"][...]
number_of_groups = len(offsets)
group_sizes = calculate_group_sizes_array(offsets, ids_in_halos.size)
group_array = create_group_array(group_sizes)
# We can now load the particle data from the snapshot.
with h5py.File("%s.hdf5" % ARGS["snapshot"], "r") as handle:
gas_particles = handle["PartType0"]
particle_ids = gas_particles["ParticleIDs"][...]
# Requires numpy 1.15 or greater.
_, particles_in_halos_mask, group_array_mask = np.intersect1d(
particle_ids,
ids_in_halos,
assume_unique=True,
return_indices=True,
)
# We also need to re-index the group array to cut out DM particles
group_array = group_array[group_array_mask]
# Kill the spare
del particle_ids
# Now we can only read the properties that we require from the snapshot!
temperatures = np.take(gas_particles["InternalEnergy"], particles_in_halos_mask)
# This 1e10 should probably be explained somewhere...
temperatures *= 1e10 * (gamma - 1) * mu * mH / kB
densities = np.take(gas_particles["Density"], particles_in_halos_mask)
# Just a quick check to make sure nothing's gone wrong.
assert len(group_array) == len(temperatures)
# Now we can loop through all the particles and find out the mean temperature and
# density in each halo.
particles_in_group = np.zeros(number_of_groups, dtype=int)
temp_in_group = np.zeros(number_of_groups, dtype=float)
dens_in_group = np.zeros(number_of_groups, dtype=float)
for group, T, rho in zip(group_array, temperatures, densities):
particles_in_group[group] += 1
temp_in_group[group] += T
dens_in_group[group] += rho
# First get a mask to ensure no runtime warnings
mask = particles_in_group != 0
# Normalize
temp_in_group[mask] /= particles_in_group[mask]
dens_in_group[mask] /= particles_in_group[mask]
# Now we can load the data according to the halo finder to compare with.
with h5py.File("%s.properties" % ARGS["velociraptor"], "r") as handle:
halo_masses = handle["Mass_200crit"][...]
halo_temperatures = virial_temp(mu, halo_masses * 1e10, h=h, a=a)
# Finally, the plotting!