Merge branch 'add-anarchy-pu' into 'master'

Adds ANARCHY-PU Scheme

See merge request !723

configure.ac

@@ -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])
+   ;;
 
 
    *)

doc/RTD/source/SubgridModels/EAGLE/index.rst

@@ -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:
 

examples/SantaBarbara/make_plots.sh

+python3 makeImage.py santabarbara_0153.hdf5 0 twilight white
+python3 plotSolution.py 153 halo
+python3 plotTempEvolution.py
+python3 rhoTHaloComparison.py

examples/SantaBarbara/plotSmoothingLength.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"])

examples/SantaBarbara/plotSolution.py

@@ -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]))
 

examples/SantaBarbara/plotTempEvolution.py

+################################################################################
+# 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)
+

examples/SantaBarbara/rhoTHaloComparison.py

+"""
+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)
+