examples/ChemistryTests/MetalAdvectionTestGEAR/makeIC.py

+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2024 Yolan Uyttenhove (yolan.uyttenhove@ugent.be)
+#
+# 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/>.
+#
+##############################################################################
+import math
+
+# ---------------------------------------------------------------------
+# Test the diffusion/advection of metals by creating regions with/without
+# initial metallicities. Run with EAGLE chemistry.
+# ---------------------------------------------------------------------
+
+import numpy as np
+import h5py
+
+GAMMA = 5 / 3
+RHO = 1
+P = 1
+VELOCITY = 2.5
+# Set ELEMENT_COUNT to 2 for AGORA runs, or for GEAR compile swift with `--with-chemistry=GEAR_X` with X equal to
+# ELEMENT_COUNT.
+ELEMENT_COUNT = 4
+
+outputfilename = "advect_metals.hdf5"
+
+
+def get_mask(boxsize, pos, element_idx):
+    x = boxsize[0]
+    right_half_mask = pos[:, 0] > x / 3
+
+    if element_idx == ELEMENT_COUNT - 1:
+        return right_half_mask
+    else:
+        periods = 2 * (element_idx + 2)
+        periods_mask = (pos[:, 0] * periods // x) % 2 == 0
+        return right_half_mask & periods_mask
+
+
+def get_abundance(element_idx):
+    if element_idx == ELEMENT_COUNT - 1:
+        return 0.2
+    else:
+        return 0.1 * 0.5 ** element_idx
+
+
+def get_element_abundances_metallicity(pos, boxsize):
+    element_abundances = []
+    for i in range(ELEMENT_COUNT):
+        element_abundances.append(
+            np.where(get_mask(boxsize, pos, i), get_abundance(i), 0)
+        )
+
+    return np.stack(element_abundances, axis=1)
+
+
+if __name__ == "__main__":
+    glass = h5py.File("glassPlane_64.hdf5", "r")
+    parts = glass["PartType0"]
+    pos = parts["Coordinates"][:]
+    pos = np.concatenate([pos, pos + np.array([1, 0, 0])])
+    h = parts["SmoothingLength"][:]
+    h = np.concatenate([h, h])
+    glass.close()
+
+    # Set up metadata
+    boxsize = np.array([2.0, 1.0, 1.0])
+    n_part = len(h)
+    ids = np.arange(n_part) + 1
+
+    # Setup other particle quantities
+    rho = RHO * np.ones_like(h)
+    rho[pos[:, 1] < 0.5 * boxsize[1]] *= 0.5
+    masses = rho * np.prod(boxsize) / n_part
+    velocities = np.zeros((n_part, 3))
+    velocities[:, :] = 0.5 * math.sqrt(2) * VELOCITY * np.array([1.0, 1.0, 0.0])
+    internal_energy = P / (rho * (GAMMA - 1))
+
+    # Setup metallicities
+    metallicities = get_element_abundances_metallicity(pos, boxsize)
+
+    # Now open the file and write the data.
+    file = h5py.File(outputfilename, "w")
+
+    # Header
+    grp = file.create_group("/Header")
+    grp.attrs["BoxSize"] = boxsize
+    grp.attrs["NumPart_Total"] = [n_part, 0, 0, 0, 0, 0]
+    grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+    grp.attrs["NumPart_ThisFile"] = [n_part, 0, 0, 0, 0, 0]
+    grp.attrs["Time"] = 0.0
+    grp.attrs["NumFilesPerSnapshot"] = 1
+    grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    grp.attrs["Flag_Entropy_ICs"] = 0
+    grp.attrs["Dimension"] = 2
+
+    # Units
+    grp = file.create_group("/Units")
+    grp.attrs["Unit length in cgs (U_L)"] = 1.0
+    grp.attrs["Unit mass in cgs (U_M)"] = 1.0
+    grp.attrs["Unit time in cgs (U_t)"] = 1.0
+    grp.attrs["Unit current in cgs (U_I)"] = 1.0
+    grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
+
+    # Particle group
+    grp = file.create_group("/PartType0")
+    grp.create_dataset("Coordinates", data=pos, dtype="d")
+    grp.create_dataset("Velocities", data=velocities, dtype="f")
+    grp.create_dataset("Masses", data=masses, dtype="f")
+    grp.create_dataset("SmoothingLength", data=h, dtype="f")
+    grp.create_dataset("InternalEnergy", data=internal_energy, dtype="f")
+    grp.create_dataset("ParticleIDs", data=ids, dtype="L")
+    grp.create_dataset("MetalMassFraction", data=metallicities, dtype="f")
+
+    file.close()
+
+    # close up, and we're done!
+    file.close()
+
+    if ELEMENT_COUNT == 2:
+        print(
+            f"Make sure swift was compiled with `--with-chemistry=GEAR_2` or `--with-chemistry=AGORA`"
+        )
+    else:
+        print(
+            f"Make sure swift was compiled with `--with-chemistry=GEAR_{ELEMENT_COUNT}`"
+        )

examples/ChemistryTests/MetalAdvectionTestGEAR/plotAdvectedMetals.py

+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2024 Yolan Uyttenhove (yolan.uyttenhove@ugent.be)
+#
+# 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/>.
+#
+##############################################################################
+
+import math
+
+import matplotlib.pyplot as plt
+from matplotlib.cm import ScalarMappable
+from matplotlib.colors import SymLogNorm, Normalize
+import numpy as np
+import unyt
+
+import swiftsimio
+from swiftsimio.visualisation import project_gas
+from pathlib import Path
+
+from makeIC import VELOCITY, ELEMENT_COUNT
+
+
+def plot_single(ax_mass, ax_fraction, name, title, data, mass_map, kwargs_inner):
+    mass_weighted_map = project_gas(data, project=name, **kwargs_inner["projection"])
+    ax_mass.imshow(mass_weighted_map.in_cgs().value.T, **kwargs_inner["imshow_mass"])
+    ax_mass.set_title(title)
+    ax_fraction.imshow(
+        (mass_weighted_map / mass_map).value.T, **kwargs_inner["imshow_fraction"]
+    )
+    ax_fraction.set_title(title)
+    ax_mass.axis("off")
+    ax_fraction.axis("off")
+
+
+def plot_all(fname, savename):
+    data = swiftsimio.load(fname)
+
+    # Shift coordinates to starting position
+    velocity = 0.5 * math.sqrt(2) * VELOCITY * np.array([1, 1, 0]) * unyt.cm / unyt.s
+    data.gas.coordinates -= data.metadata.time * velocity
+
+    # Add mass weighted element mass fractions to gas dataset
+    masses = data.gas.masses
+    element_mass_fractions = data.gas.metal_mass_fractions
+    columns = getattr(element_mass_fractions, "named_columns", None)
+    for i in range(ELEMENT_COUNT):
+        if columns is None:
+            data.gas.__setattr__(
+                f"element_mass_sp{i}", element_mass_fractions[:, i] * masses
+            )
+        else:
+            data.gas.__setattr__(
+                f"element_mass_sp{i}",
+                getattr(element_mass_fractions, columns[i]) * masses,
+            )
+
+    # Create necessary figures and axes
+    fig = plt.figure(layout="constrained", figsize=(8, 2 * ELEMENT_COUNT + 1))
+    fig.suptitle(
+        f"Profiles shifted to starting position after t={data.metadata.time:.2f}",
+        fontsize=14,
+    )
+    fig_ratios, fig_masses = fig.subfigures(1, 2)
+    fig_ratios.suptitle("Mass ratio of elements")
+    fig_masses.suptitle("Surface density in elements")
+    axes_ratios = fig_ratios.subplots(ELEMENT_COUNT, 1, sharex=True, sharey=True)
+    axes_masses = fig_masses.subplots(ELEMENT_COUNT, 1, sharex=True, sharey=True)
+
+    # parameters for swiftsimio projections
+    projection_kwargs = {
+        "region": np.array([0, 2, 0, 1, 0, 1]) * unyt.cm,
+        "resolution": 500,
+        "parallel": True,
+    }
+    # Parameters for imshow
+    if ELEMENT_COUNT > 5:
+        thresh = 10 ** math.floor(math.log10(0.5 ** (ELEMENT_COUNT - 2)))
+        norm_ratios = SymLogNorm(vmin=0, vmax=0.21, linthresh=thresh, base=10)
+        norm_masses = SymLogNorm(vmin=0, vmax=0.25, linthresh=thresh, base=10)
+    else:
+        norm_ratios = Normalize(vmin=0, vmax=0.21)
+        norm_masses = Normalize(vmin=0, vmax=0.25)
+    imshow_fraction_kwargs = dict(norm=norm_ratios, cmap="rainbow")
+    imshow_mass_kwargs = dict(norm=norm_masses, cmap="turbo")
+
+    # Plot the quantities
+    mass_map = project_gas(data, project="masses", **projection_kwargs)
+
+    # Parameters for plotting:
+    plotting_kwargs = dict(
+        data=data,
+        mass_map=mass_map,
+        kwargs_inner=dict(
+            projection=projection_kwargs,
+            imshow_mass=imshow_mass_kwargs,
+            imshow_fraction=imshow_fraction_kwargs,
+        ),
+    )
+
+    if columns is None:
+        columns = [f"Species {i + 1}" for i in range(ELEMENT_COUNT)]
+    for i in range(ELEMENT_COUNT):
+        plot_single(
+            axes_masses[i],
+            axes_ratios[i],
+            f"element_mass_sp{i}",
+            columns[i],
+            **plotting_kwargs,
+        )
+
+    # Add Colorbars
+    cb_masses = fig_masses.colorbar(
+        ScalarMappable(**imshow_mass_kwargs),
+        orientation="horizontal",
+        shrink=0.75,
+        pad=0.01,
+        ax=axes_masses,
+    )
+    cb_ratios = fig_ratios.colorbar(
+        ScalarMappable(**imshow_fraction_kwargs),
+        orientation="horizontal",
+        shrink=0.75,
+        pad=0.01,
+        ax=axes_ratios,
+    )
+    cb_masses.ax.set_xlabel("Surface density (g/cm^2)")
+    cb_ratios.ax.set_xlabel("Mass ratio")
+
+    # Save output
+    fig.savefig(savename, dpi=300)
+
+
+if __name__ == "__main__":
+    cwd = Path(__file__).parent
+    print("Plotting metals for output_0000.hdf5...")
+    plot_all(cwd / "output_0000.hdf5", savename=cwd / "metal_advection_0000.png")
+    print("Plotting metals for output_0001.hdf5...")
+    plot_all(cwd / "output_0001.hdf5", savename=cwd / "metal_advection_0001.png")
+    print("Done!")

examples/ChemistryTests/MetalAdvectionTestGEAR/run.sh

+#!/bin/bash
+
+# make run.sh fail if a subcommand fails
+set -e
+set -o pipefail
+
+if [ ! -e glassPlane_64.hdf5 ]
+then
+    echo "Fetching initial glass file ..."
+    ./getGlass.sh
+fi
+
+# Always generate ICs in case ELEMENT_COUNT has been changed
+echo "Generating ICs"
+python3 makeIC.py
+
+# Run SWIFT (must be compiled with --with-chemistry=GEAR_X or --with-chemistry=AGORA)
+../../../swift \
+    --hydro --threads=4 advect_metals.yml 2>&1 | tee output.log
+
+python3 runSanityChecksAdvectedMetals.py
+python3 plotAdvectedMetals.py

examples/ChemistryTests/MetalAdvectionTestGEAR/runSanityChecksAdvectedMetals.py

+#!/usr/bin/env python3
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2024 Yolan Uyttenhove (yolan.uyttenhove@ugent.be)
+#
+# 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/>.
+#
+##############################################################################
+
+from pathlib import Path
+
+import numpy as np
+import swiftsimio
+
+from makeIC import ELEMENT_COUNT
+
+
+def test(name):
+    def test_decorator(test_func):
+        def test_runner(*args, **kwargs):
+            try:
+                test_func(*args, **kwargs)
+            except Exception as e:
+                print(f"\tTest {name} failed!")
+                raise e
+            else:
+                print(f"\tTest {name} \u2705")
+
+        return test_runner
+
+    return test_decorator
+
+
+def approx_equal(a, b, threshold=0.001):
+    return 2 * abs(a - b) / (abs(a) + abs(b)) < threshold
+
+
+@test("total metal mass conservation")
+def test_total_metal_mass_conservation(data_start, data_end):
+    def metal_mass(data):
+        columns = getattr(data.gas.metal_mass_fractions, "named_columns", None)
+        if columns is not None:
+            return sum(
+                data.gas.masses * getattr(data.gas.metal_mass_fractions, c)
+                for c in columns
+            )
+        else:
+            return data.gas.masses.reshape(-1, 1) * data.gas.metal_mass_fractions
+
+    assert approx_equal(np.sum(metal_mass(data_start)), np.sum(metal_mass(data_end)))
+
+
+def element_mass(data, element_idx):
+    columns = getattr(data.gas.metal_mass_fractions, "named_columns", None)
+    if columns is not None:
+        return data.gas.masses * getattr(
+            data.gas.metal_mass_fractions, columns[element_idx]
+        )
+    else:
+        return data.gas.masses * data.gas.metal_mass_fractions[:, element_idx]
+
+
+@test("element-wise mass conservation")
+def test_element_wise_mass_conservation(data_start, data_end):
+    for i in range(ELEMENT_COUNT):
+        assert approx_equal(
+            np.sum(element_mass(data_start, i)), np.sum(element_mass(data_end, i))
+        )
+
+
+if __name__ == "__main__":
+    print("Running sanity checks...")
+
+    cwd = Path(__file__).parent
+    start = swiftsimio.load(cwd / "output_0000.hdf5")
+    end = swiftsimio.load(cwd / "output_0001.hdf5")
+
+    test_total_metal_mass_conservation(start, end)
+    test_element_wise_mass_conservation(start, end)
+
+    print("Done!")

examples/Cooling/ConstantCosmoTempEvolution/const_cosmo_temp_evol.yml

@@ -11,7 +11,7 @@ Cosmology:
   h:              0.6777        # Reduced Hubble constant
   a_begin:        0.0099        # Initial scale-factor of the simulation (z = 100.0)
   a_end:          1.0           # Final scale factor of the simulation
-  Omega_m:        0.307         # Matter density parameter
+  Omega_cdm:      0.2587481     # Cold Dark Matter density parameter
   Omega_lambda:   0.693         # Dark-energy density parameter
   Omega_b:        0.0482519     # Baryon density parameter
 

examples/Cooling/ConstantCosmoTempEvolution/getGlass.sh

 #!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/gravity_glassCube_32.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/gravity_glassCube_32.hdf5

examples/Cooling/ConstantCosmoTempEvolution/makeIC.py

@@ -51,7 +51,7 @@ glass.close()
 # Calculate mean baryon density today from comological parameters
 H_0 = 100.0 * hubble_param * unyt.km / unyt.s / unyt.Mpc
 rho_crit_0 = 3.0 * H_0 ** 2 / (8.0 * np.pi * unyt.G)
-rho_bar_0 = Omega_bar * rho_crit_0 
+rho_bar_0 = Omega_bar * rho_crit_0
 
 # From this, we calculate the mass of the gas particles
 gas_particle_mass = rho_bar_0 * boxsize ** 3 / (n_p_1D ** 3)

examples/Cooling/ConstantCosmoTempEvolution/plot_thermal_history.py

@@ -84,8 +84,8 @@ for index, data in enumerate(snap_data):
 # first calculate rho_bar_0 from snapshot metadata
 data = snap_data[0]
 cosmology = data.metadata.cosmology
-H0 = cosmology["H0 [internal units]"] / (data.units.time)
-Omega_bar = cosmology["Omega_b"]
+H0 = unyt.unyt_quantity.from_astropy(cosmology.H0)
+Omega_bar = unyt.unyt_quantity(cosmology.Ob0)
 
 ### now calculate rho_bar_0 and divide through
 rho_bar_0 = 3.0 * H0 ** 2 / (8.0 * np.pi * unyt.G) * Omega_bar

examples/Cooling/ConstantCosmoTempEvolution/run.sh

@@ -20,7 +20,7 @@ then
 fi
 
 # Run SWIFT
-../../swift --hydro --cosmology --cooling --threads=4 const_cosmo_temp_evol.yml 2>&1 | tee output.log
+../../../swift --hydro --cosmology --cooling --threads=4 const_cosmo_temp_evol.yml 2>&1 | tee output.log
 
 # Plot the result
 python3 plot_thermal_history.py cooling_box

examples/Cooling/CoolingBox/README

+Runs a uniform box exposed to constant cooling rates computed with 
+the Grackle library.
+
+To run with Grackle mode 0, configure swift with:
+  ./configure --with-cooling=grackle_0 --with-grackle=$GRACKLE_ROOT
+
+To run with Grackle mode 1, configure swift with:
+  ./configure --with-cooling=grackle_1 --with-grackle=$GRACKLE_ROOT
+  
+To run with Grackle mode 2, configure swift with:
+  ./configure --with-cooling=grackle_2 --with-grackle=$GRACKLE_ROOT
+  
+To run with Grackle mode 3, configure swift with:
+  ./configure --with-cooling=grackle_3 --with-grackle=$GRACKLE_ROOT    

examples/Cooling/CoolingBox/coolingBox.yml

@@ -23,6 +23,9 @@ Snapshots:
 Statistics:
   delta_time:          1e-3 # Time between statistics output
 
+Scheduler:
+  tasks_per_cell: 64
+
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
@@ -52,6 +55,7 @@ GrackleCooling:
   thermal_time_myr: 5                          # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
   self_shielding_method: 0                    # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
   self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
+  maximal_density_Hpcm3: -1                    # Maximal density (in hydrogen atoms/cm^3) for cooling. Higher densities are floored to this value to ensure grackle works properly when interpolating beyond the cloudy_table maximal density. A value < 0 deactivates this parameter.
 
 GEARChemistry:
   initial_metallicity: 0.01295
@@ -78,3 +82,13 @@ EAGLEChemistry:             # Solar abundances
 
 GEARPressureFloor:
   jeans_factor: 0.       # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+
+PS2020Cooling:
+  dir_name:                ./UV_dust1_CR1_G1_shield1.hdf5 # Location of the cooling tables
+  H_reion_z:               7.5               # Redshift of Hydrogen re-ionization (Planck 2018)
+  H_reion_eV_p_H:          2.0
+  He_reion_z_centre:       3.5               # Redshift of the centre of the Helium re-ionization Gaussian
+  He_reion_z_sigma:        0.5               # Spread in redshift of the  Helium re-ionization Gaussian
+  He_reion_eV_p_H:         2.0               # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
+  delta_logTEOS_subgrid_properties: 0.3      # delta log T above the EOS below which the subgrid properties use Teq assumption
+  rapid_cooling_threshold:          0.333333 # Switch to rapid cooling regime for dt / t_cool above this threshold.

examples/Cooling/CoolingBox/getGlass.sh

 #!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5

examples/Cooling/CoolingBox/makeIC.py

 ###############################################################################
 # This file is part of SWIFT.
 # Copyright (c) 2016 Stefan Arridge (stefan.arridge@durhama.ac.uk)
-#                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#                    Matthieu Schaller (schaller@strw.leidenuniv.nl)
 #
 # 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
@@ -28,14 +28,14 @@ periodic = 1  # 1 For periodic box
 boxSize = 1  # 1 kiloparsec
 rho = 3.2e3  # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
 T = 4e3  # Initial Temperature
-gamma = 5./3.  # Gas adiabatic index
+gamma = 5.0 / 3.0  # Gas adiabatic index
 fileName = "coolingBox.hdf5"
 # ---------------------------------------------------
 
 # defines some constants
 # need to be changed in plotTemperature.py too
 h_frac = 0.76
-mu = 4. / (1. + 3. * h_frac)
+mu = 4.0 / (1.0 + 3.0 * h_frac)
 
 m_h_cgs = 1.67e-24
 k_b_cgs = 1.38e-16
@@ -53,12 +53,12 @@ h = glass["/PartType0/SmoothingLength"][:] * boxSize
 
 # Compute basic properties
 numPart = np.size(pos) // 3
-mass = boxSize**3 * rho / numPart
-internalEnergy = k_b_cgs * T * mu / ((gamma - 1.) * m_h_cgs)
-internalEnergy *= (unit_time / unit_length)**2
+mass = boxSize ** 3 * rho / numPart
+internalEnergy = k_b_cgs * T * mu / ((gamma - 1.0) * m_h_cgs)
+internalEnergy *= (unit_time / unit_length) ** 2
 
 # File
-f = h5py.File(fileName, 'w')
+f = h5py.File(fileName, "w")
 
 # Header
 grp = f.create_group("/Header")
@@ -80,33 +80,33 @@ grp = f.create_group("/Units")
 grp.attrs["Unit length in cgs (U_L)"] = unit_length
 grp.attrs["Unit mass in cgs (U_M)"] = unit_mass
 grp.attrs["Unit time in cgs (U_t)"] = unit_time
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
 
 # Particle group
 grp = f.create_group("/PartType0")
 
 v = np.zeros((numPart, 3))
-ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ds = grp.create_dataset("Velocities", (numPart, 3), "f")
 ds[()] = v
 
 m = np.full((numPart, 1), mass)
-ds = grp.create_dataset('Masses', (numPart, 1), 'f')
+ds = grp.create_dataset("Masses", (numPart, 1), "f")
 ds[()] = m
 
 h = np.reshape(h, (numPart, 1))
-ds = grp.create_dataset('SmoothingLength', (numPart, 1), 'f')
+ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
 ds[()] = h
 
 u = np.full((numPart, 1), internalEnergy)
-ds = grp.create_dataset('InternalEnergy', (numPart, 1), 'f')
+ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
 ds[()] = u
 
 ids = np.reshape(ids, (numPart, 1))
-ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
+ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
 ds[()] = ids
 
-ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
 ds[()] = pos
 
 f.close()

examples/Cooling/CoolingBox/plotEnergy.py

-from h5py import File
-import numpy as np
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2022 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+#
+# 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/>.
+#
+##############################################################################
 import matplotlib
-from glob import glob
+
 matplotlib.use("Agg")
 import matplotlib.pyplot as plt
+import numpy as np
+from glob import glob
+import h5py
 
-# Plot parameters
-params = {
-    'axes.labelsize': 10,
-    'axes.titlesize': 10,
-    'font.size': 12,
-    'legend.fontsize': 12,
-    'xtick.labelsize': 10,
-    'ytick.labelsize': 10,
-    'text.usetex': True,
-    'figure.figsize': (5, 5),
-    'figure.subplot.left': 0.145,
-    'figure.subplot.right': 0.99,
-    'figure.subplot.bottom': 0.11,
-    'figure.subplot.top': 0.99,
-    'figure.subplot.wspace': 0.15,
-    'figure.subplot.hspace': 0.12,
-    'lines.markersize': 6,
-    'lines.linewidth': 3.,
-}
-plt.rcParams.update(params)
-
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 # Some constants in cgs units
 k_b_cgs = 1.38e-16  # boltzmann
 m_h_cgs = 1.67e-24  # proton mass
 
-
 # File containing the total energy
-stats_filename = "./energy.txt"
+stats_filename = "./statistics.txt"
 
 # First snapshot
 snap_filename = "coolingBox_0000.hdf5"
 
 # Read the initial state of the gas
-f = File(snap_filename, 'r')
+f = h5py.File(snap_filename, "r")
 
 # Read the units parameters from the snapshot
 units = f["InternalCodeUnits"]
@@ -53,18 +51,19 @@ gamma = float(f["HydroScheme"].attrs["Adiabatic index"])
 
 def energyUnits(u):
     """ Compute the temperature from the internal energy. """
-    u *= (unit_length / unit_time)**2
+    u *= (unit_length / unit_time) ** 2
     return u * m_h_cgs / k_b_cgs
 
 
 # Read energy and time arrays
 array = np.genfromtxt(stats_filename, skip_header=1)
-time = array[:, 0] * unit_time
-total_mass = array[:, 1]
-total_energy = array[:, 2]
-kinetic_energy = array[:, 3]
-internal_energy = array[:, 4]
-radiated_energy = array[:, 8]
+time = array[:, 1] * unit_time
+total_mass = array[:, 4]
+kinetic_energy = array[:, 13]
+internal_energy = array[:, 14]
+potential_energy = array[:, 15]
+radiated_energy = array[:, 16]
+total_energy = kinetic_energy + internal_energy + potential_energy
 initial_energy = total_energy[0]
 
 # Conversions to cgs
@@ -86,7 +85,7 @@ N = len(files)
 temp_snap = np.zeros(N)
 time_snap_cgs = np.zeros(N)
 for i in range(N):
-    snap = File(files[i], 'r')
+    snap = h5py.File(files[i], "r")
     u = snap["/PartType0/InternalEnergies"][:] * snap["/PartType0/Masses"][:]
     u = sum(u) / total_mass[0]
     temp_snap[i] = energyUnits(u)
@@ -95,16 +94,23 @@ for i in range(N):
 
 plt.figure()
 
-Myr_in_yr = 3.15e13
-plt.plot(time, total_energy_cgs, 'r-', lw=1.6, label="Gas total energy")
-plt.plot(time_snap_cgs, temp_snap, 'rD', ms=3)
-plt.plot(time, radiated_energy_cgs, 'g-', lw=1.6, label="Radiated energy")
-plt.plot(time, total_energy_cgs + radiated_energy_cgs, 'b-',
-         lw=0.6, label="Gas total + radiated")
-
-plt.legend(loc="right", fontsize=8, frameon=False,
-           handlelength=3, ncol=1)
+Myr_in_s = 1e6 * 365.25 * 24.0 * 60.0 * 60.0
+plt.plot(time / Myr_in_s, total_energy_cgs, "r-", lw=1.6, label="Gas total energy")
+# statistics and snapshots may not be at same timestep and frequency
+plt.plot(time_snap_cgs / Myr_in_s, temp_snap, "rD", ms=3)
+plt.plot(time / Myr_in_s, radiated_energy_cgs, "g-", lw=1.6, label="Radiated energy")
+plt.plot(
+    time / Myr_in_s,
+    total_energy_cgs + radiated_energy_cgs,
+    "b-",
+    lw=0.6,
+    label="Gas total + radiated",
+)
+
+plt.legend(loc="right", fontsize=8, frameon=False, handlelength=3, ncol=1)
 plt.xlabel("${\\rm{Time~[Myr]}}$", labelpad=0)
 plt.ylabel("${\\rm{Internal ~Energy ~(u ~m_H / k_B) ~[K]}}$")
 
+plt.tight_layout()
+
 plt.savefig("energy.png", dpi=200)

examples/Cooling/CoolingBox/run.sh

@@ -10,18 +10,18 @@ fi
 if [ ! -e coolingBox.hdf5 ]
 then
     echo "Generating initial conditions for the cooling box example..."
-    python makeIC.py
+    python3 makeIC.py
 fi
 
 # Get the Grackle cooling table
 if [ ! -e CloudyData_UVB=HM2012.h5 ]
 then
     echo "Fetching the Cloudy tables required by Grackle..."
-    ../getCoolingTable.sh
+    ../getGrackleCoolingTable.sh
 fi
 
 # Run SWIFT
-../../swift --hydro --cooling --threads=4 -n 1000 coolingBox.yml
+../../../swift --hydro --cooling --threads=4 -n 1000 coolingBox.yml
 
 # Check energy conservation and cooling rate
-python plotEnergy.py
+python3 plotEnergy.py

examples/Cooling/CoolingHalo/density_profile.py

@@ -5,23 +5,23 @@ import sys
 
 n_snaps = 11
 
-#for the plotting
-#n_radial_bins = int(sys.argv[1])
+# for the plotting
+# n_radial_bins = int(sys.argv[1])
 
-#some constants
-OMEGA = 0.3 # Cosmological matter fraction at z = 0
+# some constants
+OMEGA = 0.3  # Cosmological matter fraction at z = 0
 PARSEC_IN_CGS = 3.0856776e18
 KM_PER_SEC_IN_CGS = 1.0e5
 CONST_G_CGS = 6.672e-8
-h = 0.67777 # hubble parameter
-gamma = 5./3.
+h = 0.67777  # hubble parameter
+gamma = 5.0 / 3.0
 eta = 1.2349
-H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+H_0_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
 
-#read some header/parameter information from the first snapshot
+# read some header/parameter information from the first snapshot
 
 filename = "Hydrostatic_0000.hdf5"
-f = h5.File(filename,'r')
+f = h5.File(filename, "r")
 params = f["Parameters"]
 unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
@@ -33,69 +33,71 @@ header = f["Header"]
 N = header.attrs["NumPart_Total"][0]
 box_centre = np.array(header.attrs["BoxSize"])
 
-#calculate r_vir and M_vir from v_c
-r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
-M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+# calculate r_vir and M_vir from v_c
+r_vir_cgs = v_c_cgs / (10.0 * H_0_cgs * np.sqrt(OMEGA))
+M_vir_cgs = r_vir_cgs * v_c_cgs ** 2 / CONST_G_CGS
 
 for i in range(n_snaps):
 
-    filename = "Hydrostatic_%04d.hdf5" %i
-    f = h5.File(filename,'r')
+    filename = "Hydrostatic_%04d.hdf5" % i
+    f = h5.File(filename, "r")
     coords_dset = f["PartType0/Coordinates"]
     coords = np.array(coords_dset)
-#translate coords by centre of box
+    # translate coords by centre of box
     header = f["Header"]
     snap_time = header.attrs["Time"]
     snap_time_cgs = snap_time * unit_time_cgs
-    coords[:,0] -= box_centre[0]/2.
-    coords[:,1] -= box_centre[1]/2.
-    coords[:,2] -= box_centre[2]/2.
-    radius = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2)
-    radius_cgs = radius*unit_length_cgs
+    coords[:, 0] -= box_centre[0] / 2.0
+    coords[:, 1] -= box_centre[1] / 2.0
+    coords[:, 2] -= box_centre[2] / 2.0
+    radius = np.sqrt(coords[:, 0] ** 2 + coords[:, 1] ** 2 + coords[:, 2] ** 2)
+    radius_cgs = radius * unit_length_cgs
     radius_over_virial_radius = radius_cgs / r_vir_cgs
 
     r = radius_over_virial_radius
 
     # bin_width = 1./n_radial_bins
-#     hist = np.histogram(r,bins = n_radial_bins)[0] # number of particles in each bin
+    #     hist = np.histogram(r,bins = n_radial_bins)[0] # number of particles in each bin
 
-# #find the mass in each radial bin
+    # #find the mass in each radial bin
 
-#     mass_dset = f["PartType0/Masses"]
-# #mass of each particles should be equal
-#     part_mass = np.array(mass_dset)[0]
-#     part_mass_cgs = part_mass * unit_mass_cgs
-#     part_mass_over_virial_mass = part_mass_cgs / M_vir_cgs 
+    #     mass_dset = f["PartType0/Masses"]
+    # #mass of each particles should be equal
+    #     part_mass = np.array(mass_dset)[0]
+    #     part_mass_cgs = part_mass * unit_mass_cgs
+    #     part_mass_over_virial_mass = part_mass_cgs / M_vir_cgs
 
-#     mass_hist = hist * part_mass_over_virial_mass
-#     radial_bin_mids = np.linspace(bin_width/2.,1 - bin_width/2.,n_radial_bins)
-# #volume in each radial bin
-#     volume = 4.*np.pi * radial_bin_mids**2 * bin_width
+    #     mass_hist = hist * part_mass_over_virial_mass
+    #     radial_bin_mids = np.linspace(bin_width/2.,1 - bin_width/2.,n_radial_bins)
+    # #volume in each radial bin
+    #     volume = 4.*np.pi * radial_bin_mids**2 * bin_width
 
-# #now divide hist by the volume so we have a density in each bin
+    # #now divide hist by the volume so we have a density in each bin
 
-#     density = mass_hist / volume
+    #     density = mass_hist / volume
 
     # read the densities
 
     density_dset = f["PartType0/Density"]
     density = np.array(density_dset)
-    density_cgs = density * unit_mass_cgs / unit_length_cgs**3
-    rho = density_cgs * r_vir_cgs**3 / M_vir_cgs
+    density_cgs = density * unit_mass_cgs / unit_length_cgs ** 3
+    rho = density_cgs * r_vir_cgs ** 3 / M_vir_cgs
 
-    t = np.linspace(0.01,2.0,1000)
-    rho_analytic = t**(-2)/(4.*np.pi)
+    t = np.linspace(0.01, 2.0, 1000)
+    rho_analytic = t ** (-2) / (4.0 * np.pi)
 
-    plt.plot(r,rho,'x',label = "Numerical solution")
-    plt.plot(t,rho_analytic,label = "Analytic Solution")
-    plt.legend(loc = "upper right")
+    plt.plot(r, rho, "x", label="Numerical solution")
+    plt.plot(t, rho_analytic, label="Analytic Solution")
+    plt.legend(loc="upper right")
     plt.xlabel(r"$r / r_{vir}$")
     plt.ylabel(r"$\rho / (M_{vir} / r_{vir}^3)$")
-    plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    #plt.ylim((0.1,40))
-    plt.xscale('log')
-    plt.yscale('log')
-    plot_filename = "density_profile_%03d.png" %i
-    plt.savefig(plot_filename,format = "png")
+    plt.title(
+        r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$"
+        % (snap_time_cgs, N, v_c)
+    )
+    # plt.ylim((0.1,40))
+    plt.xscale("log")
+    plt.yscale("log")
+    plot_filename = "density_profile_%03d.png" % i
+    plt.savefig(plot_filename, format="png")
     plt.close()
-

examples/Cooling/CoolingHalo/internal_energy_profile.py

@@ -3,43 +3,46 @@ import h5py as h5
 import matplotlib.pyplot as plt
 import sys
 
-def do_binning(x,y,x_bin_edges):
 
-    #x and y are arrays, where y = f(x)
-    #returns number of elements of x in each bin, and the total of the y elements corresponding to those x values
+def do_binning(x, y, x_bin_edges):
+
+    # x and y are arrays, where y = f(x)
+    # returns number of elements of x in each bin, and the total of the y elements corresponding to those x values
 
     n_bins = x_bin_edges.size - 1
     count = np.zeros(n_bins)
     y_totals = np.zeros(n_bins)
-    
+
     for i in range(n_bins):
-        ind = np.intersect1d(np.where(x > bin_edges[i])[0],np.where(x <= bin_edges[i+1])[0])
+        ind = np.intersect1d(
+            np.where(x > bin_edges[i])[0], np.where(x <= bin_edges[i + 1])[0]
+        )
         count[i] = ind.size
         binned_y = y[ind]
         y_totals[i] = np.sum(binned_y)
 
-    return(count,y_totals)
+    return (count, y_totals)
 
 
 n_snaps = 100
 
-#for the plotting
+# for the plotting
 n_radial_bins = int(sys.argv[1])
 
-#some constants
-OMEGA = 0.3 # Cosmological matter fraction at z = 0
+# some constants
+OMEGA = 0.3  # Cosmological matter fraction at z = 0
 PARSEC_IN_CGS = 3.0856776e18
 KM_PER_SEC_IN_CGS = 1.0e5
 CONST_G_CGS = 6.672e-8
-h = 0.67777 # hubble parameter
-gamma = 5./3.
+h = 0.67777  # hubble parameter
+gamma = 5.0 / 3.0
 eta = 1.2349
-H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+H_0_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
 
-#read some header/parameter information from the first snapshot
+# read some header/parameter information from the first snapshot
 
 filename = "Hydrostatic_0000.hdf5"
-f = h5.File(filename,'r')
+f = h5.File(filename, "r")
 params = f["Parameters"]
 unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
@@ -51,54 +54,54 @@ header = f["Header"]
 N = header.attrs["NumPart_Total"][0]
 box_centre = np.array(header.attrs["BoxSize"])
 
-#calculate r_vir and M_vir from v_c
-r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
-M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+# calculate r_vir and M_vir from v_c
+r_vir_cgs = v_c_cgs / (10.0 * H_0_cgs * np.sqrt(OMEGA))
+M_vir_cgs = r_vir_cgs * v_c_cgs ** 2 / CONST_G_CGS
 
 for i in range(n_snaps):
 
-    filename = "Hydrostatic_%04d.hdf5" %i
-    f = h5.File(filename,'r')
+    filename = "Hydrostatic_%04d.hdf5" % i
+    f = h5.File(filename, "r")
     coords_dset = f["PartType0/Coordinates"]
     coords = np.array(coords_dset)
-#translate coords by centre of box
+    # translate coords by centre of box
     header = f["Header"]
     snap_time = header.attrs["Time"]
     snap_time_cgs = snap_time * unit_time_cgs
-    coords[:,0] -= box_centre[0]/2.
-    coords[:,1] -= box_centre[1]/2.
-    coords[:,2] -= box_centre[2]/2.
-    radius = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2)
-    radius_cgs = radius*unit_length_cgs
+    coords[:, 0] -= box_centre[0] / 2.0
+    coords[:, 1] -= box_centre[1] / 2.0
+    coords[:, 2] -= box_centre[2] / 2.0
+    radius = np.sqrt(coords[:, 0] ** 2 + coords[:, 1] ** 2 + coords[:, 2] ** 2)
+    radius_cgs = radius * unit_length_cgs
     radius_over_virial_radius = radius_cgs / r_vir_cgs
 
-#get the internal energies
+    # get the internal energies
     u_dset = f["PartType0/InternalEnergy"]
     u = np.array(u_dset)
 
-#make dimensionless
-    u /= v_c**2/(2. * (gamma - 1.))
+    # make dimensionless
+    u /= v_c ** 2 / (2.0 * (gamma - 1.0))
     r = radius_over_virial_radius
 
-    bin_edges = np.linspace(0,1,n_radial_bins + 1)
-    (hist,u_totals) = do_binning(r,u,bin_edges)
-    
-    bin_widths = 1. / n_radial_bins
-    radial_bin_mids = np.linspace(bin_widths / 2. , 1. - bin_widths / 2. , n_radial_bins) 
-    binned_u = u_totals / hist
+    bin_edges = np.linspace(0, 1, n_radial_bins + 1)
+    (hist, u_totals) = do_binning(r, u, bin_edges)
 
+    bin_widths = 1.0 / n_radial_bins
+    radial_bin_mids = np.linspace(
+        bin_widths / 2.0, 1.0 - bin_widths / 2.0, n_radial_bins
+    )
+    binned_u = u_totals / hist
 
-    plt.plot(radial_bin_mids,binned_u,'ko',label = "Numerical solution")
-    plt.plot((0,1),(1,1),label = "Analytic Solution")
-    plt.legend(loc = "lower right")
+    plt.plot(radial_bin_mids, binned_u, "ko", label="Numerical solution")
+    plt.plot((0, 1), (1, 1), label="Analytic Solution")
+    plt.legend(loc="lower right")
     plt.xlabel(r"$r / r_{vir}$")
     plt.ylabel(r"$u / (v_c^2 / (2(\gamma - 1)) $")
-    plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
-    plt.ylim((0,1))
-    plot_filename = "internal_energy_profile_%03d.png" %i
-    plt.savefig(plot_filename,format = "png")
+    plt.title(
+        r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$"
+        % (snap_time_cgs, N, v_c)
+    )
+    plt.ylim((0, 1))
+    plot_filename = "internal_energy_profile_%03d.png" % i
+    plt.savefig(plot_filename, format="png")
     plt.close()
-
-
-        
-    

examples/Cooling/CoolingHalo/makeIC.py

 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2016 Stefan Arridge (stefan.arridge@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/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Stefan Arridge (stefan.arridge@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/>.
+#
+##############################################################################
 
 import h5py
 import sys
@@ -24,71 +24,79 @@ import math
 import random
 
 # Generates N particles in a spherically symmetric distribution with density profile ~r^(-2)
-# usage: python makeIC.py 1000: generate 1000 particles
+# usage: python3 makeIC.py 1000: generate 1000 particles
 
 # Some constants
 
-OMEGA = 0.3 # Cosmological matter fraction at z = 0
+OMEGA = 0.3  # Cosmological matter fraction at z = 0
 PARSEC_IN_CGS = 3.0856776e18
 KM_PER_SEC_IN_CGS = 1.0e5
 CONST_G_CGS = 6.672e-8
-h = 0.67777 # hubble parameter
-gamma = 5./3.
+h = 0.67777  # hubble parameter
+gamma = 5.0 / 3.0
 eta = 1.2349
 
 # First set unit velocity and then the circular velocity parameter for the isothermal potential
-const_unit_velocity_in_cgs = 1.e5 #kms^-1
+const_unit_velocity_in_cgs = 1.0e5  # kms^-1
 
-v_c = 200.
+v_c = 200.0
 v_c_cgs = v_c * const_unit_velocity_in_cgs
 
 # Now we use this to get the virial mass and virial radius, which we will set to be the unit mass and radius
 
 # Find H_0, the inverse Hubble time, in cgs
 
-H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+H_0_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
 
 # From this we can find the virial radius, the radius within which the average density of the halo is
 # 200. * the mean matter density
 
-r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
+r_vir_cgs = v_c_cgs / (10.0 * H_0_cgs * np.sqrt(OMEGA))
 
 # Now get the virial mass
 
-M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+M_vir_cgs = r_vir_cgs * v_c_cgs ** 2 / CONST_G_CGS
 
 # Now set the unit length and mass
 
 const_unit_mass_in_cgs = M_vir_cgs
 const_unit_length_in_cgs = r_vir_cgs
 
-print "UnitMass_in_cgs:     ", const_unit_mass_in_cgs 
-print "UnitLength_in_cgs:   ", const_unit_length_in_cgs
-print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
-
-#derived quantities
-const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
-print "UnitTime_in_cgs:     ", const_unit_time_in_cgs
-const_G                = ((CONST_G_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
-print 'G=', const_G
+print("UnitMass_in_cgs:     ", const_unit_mass_in_cgs)
+print("UnitLength_in_cgs:   ", const_unit_length_in_cgs)
+print("UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs)
+
+# derived quantities
+const_unit_time_in_cgs = const_unit_length_in_cgs / const_unit_velocity_in_cgs
+print("UnitTime_in_cgs:     ", const_unit_time_in_cgs)
+const_G = (
+    CONST_G_CGS
+    * const_unit_mass_in_cgs
+    * const_unit_time_in_cgs
+    * const_unit_time_in_cgs
+    / (const_unit_length_in_cgs * const_unit_length_in_cgs * const_unit_length_in_cgs)
+)
+print("G=", const_G)
 
 # Parameters
-periodic= 1            # 1 For periodic box
-boxSize = 4.          
-G       = const_G 
-N       = int(sys.argv[1])  # Number of particles
+periodic = 1  # 1 For periodic box
+boxSize = 4.0
+G = const_G
+N = int(sys.argv[1])  # Number of particles
 
 # Create the file
 filename = "CoolingHalo.hdf5"
-file = h5py.File(filename, 'w')
+file = h5py.File(filename, "w")
 
-#Units
+# Units
 grp = file.create_group("/Units")
 grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
-grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs 
-grp.attrs["Unit time in cgs (U_t)"] = const_unit_length_in_cgs / const_unit_velocity_in_cgs
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs
+grp.attrs["Unit time in cgs (U_t)"] = (
+    const_unit_length_in_cgs / const_unit_velocity_in_cgs
+)
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
 
 
 # set seed for random number
@@ -97,37 +105,39 @@ np.random.seed(1234)
 
 # Positions
 # r^(-2) distribution corresponds to uniform distribution in radius
-radius = boxSize  * np.sqrt(3.) / 2.* np.random.rand(N) #the diagonal extent of the cube
-ctheta = -1. + 2 * np.random.rand(N)
-stheta = np.sqrt(1.-ctheta**2)
-phi    =  2 * math.pi * np.random.rand(N)
-coords      = np.zeros((N, 3))
-coords[:,0] = radius * stheta * np.cos(phi)
-coords[:,1] = radius * stheta * np.sin(phi)
-coords[:,2] = radius * ctheta
-
-#shift to centre of box
-coords += np.full((N,3),boxSize/2.)
-print "x range = (%f,%f)" %(np.min(coords[:,0]),np.max(coords[:,0]))
-print "y range = (%f,%f)" %(np.min(coords[:,1]),np.max(coords[:,1]))
-print "z range = (%f,%f)" %(np.min(coords[:,2]),np.max(coords[:,2]))
-
-print np.mean(coords[:,0])
-print np.mean(coords[:,1])
-print np.mean(coords[:,2])
-
-#now find the particles which are within the box
-
-x_coords = coords[:,0]
-y_coords = coords[:,1]
-z_coords = coords[:,2]
+radius = (
+    boxSize * np.sqrt(3.0) / 2.0 * np.random.rand(N)
+)  # the diagonal extent of the cube
+ctheta = -1.0 + 2 * np.random.rand(N)
+stheta = np.sqrt(1.0 - ctheta ** 2)
+phi = 2 * math.pi * np.random.rand(N)
+coords = np.zeros((N, 3))
+coords[:, 0] = radius * stheta * np.cos(phi)
+coords[:, 1] = radius * stheta * np.sin(phi)
+coords[:, 2] = radius * ctheta
+
+# shift to centre of box
+coords += np.full((N, 3), boxSize / 2.0)
+print("x range = (%f,%f)" % (np.min(coords[:, 0]), np.max(coords[:, 0])))
+print("y range = (%f,%f)" % (np.min(coords[:, 1]), np.max(coords[:, 1])))
+print("z range = (%f,%f)" % (np.min(coords[:, 2]), np.max(coords[:, 2])))
+
+print(np.mean(coords[:, 0]))
+print(np.mean(coords[:, 1]))
+print(np.mean(coords[:, 2]))
+
+# now find the particles which are within the box
+
+x_coords = coords[:, 0]
+y_coords = coords[:, 1]
+z_coords = coords[:, 2]
 
 ind = np.where(x_coords < boxSize)[0]
 x_coords = x_coords[ind]
 y_coords = y_coords[ind]
 z_coords = z_coords[ind]
 
-ind = np.where(x_coords > 0.)[0]
+ind = np.where(x_coords > 0.0)[0]
 x_coords = x_coords[ind]
 y_coords = y_coords[ind]
 z_coords = z_coords[ind]
@@ -137,7 +147,7 @@ x_coords = x_coords[ind]
 y_coords = y_coords[ind]
 z_coords = z_coords[ind]
 
-ind = np.where(y_coords > 0.)[0]
+ind = np.where(y_coords > 0.0)[0]
 x_coords = x_coords[ind]
 y_coords = y_coords[ind]
 z_coords = z_coords[ind]
@@ -147,39 +157,39 @@ x_coords = x_coords[ind]
 y_coords = y_coords[ind]
 z_coords = z_coords[ind]
 
-ind = np.where(z_coords > 0.)[0]
+ind = np.where(z_coords > 0.0)[0]
 x_coords = x_coords[ind]
 y_coords = y_coords[ind]
 z_coords = z_coords[ind]
 
-#count number of particles
+# count number of particles
 
 N = x_coords.size
 
-print "Number of particles in the box = " , N
+print("Number of particles in the box = ", N)
 
-#make the coords and radius arrays again
-coords_2 = np.zeros((N,3))
-coords_2[:,0] = x_coords
-coords_2[:,1] = y_coords
-coords_2[:,2] = z_coords
+# make the coords and radius arrays again
+coords_2 = np.zeros((N, 3))
+coords_2[:, 0] = x_coords
+coords_2[:, 1] = y_coords
+coords_2[:, 2] = z_coords
 
-radius = np.sqrt(coords_2[:,0]**2 + coords_2[:,1]**2 + coords_2[:,2]**2)
+radius = np.sqrt(coords_2[:, 0] ** 2 + coords_2[:, 1] ** 2 + coords_2[:, 2] ** 2)
 
-#test we've done it right
+# test we've done it right
 
-print "x range = (%f,%f)" %(np.min(coords_2[:,0]),np.max(coords_2[:,0]))
-print "y range = (%f,%f)" %(np.min(coords_2[:,1]),np.max(coords_2[:,1]))
-print "z range = (%f,%f)" %(np.min(coords_2[:,2]),np.max(coords_2[:,2]))
+print("x range = (%f,%f)" % (np.min(coords_2[:, 0]), np.max(coords_2[:, 0])))
+print("y range = (%f,%f)" % (np.min(coords_2[:, 1]), np.max(coords_2[:, 1])))
+print("z range = (%f,%f)" % (np.min(coords_2[:, 2]), np.max(coords_2[:, 2])))
 
-print np.mean(coords_2[:,0])
-print np.mean(coords_2[:,1])
-print np.mean(coords_2[:,2])
+print(np.mean(coords_2[:, 0]))
+print(np.mean(coords_2[:, 1]))
+print(np.mean(coords_2[:, 2]))
 
 # Header
 grp = file.create_group("/Header")
 grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] =  [N ,0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total"] = [N, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_ThisFile"] = [N, 0, 0, 0, 0, 0]
 grp.attrs["Time"] = 0.0
@@ -191,40 +201,42 @@ grp.attrs["Dimension"] = 3
 # Particle group
 grp = file.create_group("/PartType0")
 
-ds = grp.create_dataset('Coordinates', (N, 3), 'd')
+ds = grp.create_dataset("Coordinates", (N, 3), "d")
 ds[()] = coords_2
 coords_2 = np.zeros(1)
 
 # All velocities set to zero
-v = np.zeros((N,3))
-ds = grp.create_dataset('Velocities', (N, 3), 'f')
+v = np.zeros((N, 3))
+ds = grp.create_dataset("Velocities", (N, 3), "f")
 ds[()] = v
 v = np.zeros(1)
 
 # All particles of equal mass
-mass = 1. / N
-m = np.full((N,),mass)
-ds = grp.create_dataset('Masses', (N, ), 'f')
+mass = 1.0 / N
+m = np.full((N,), mass)
+ds = grp.create_dataset("Masses", (N,), "f")
 ds[()] = m
 m = np.zeros(1)
 
 # Smoothing lengths
-l = (4. * np.pi * radius**2 / N)**(1./3.) #local mean inter-particle separation
-h = np.full((N, ), eta * l)
-ds = grp.create_dataset('SmoothingLength', (N,), 'f')
+l = (4.0 * np.pi * radius ** 2 / N) ** (
+    1.0 / 3.0
+)  # local mean inter-particle separation
+h = np.full((N,), eta * l)
+ds = grp.create_dataset("SmoothingLength", (N,), "f")
 ds[()] = h
 h = np.zeros(1)
 
 # Internal energies
-u = v_c**2 / (2. * (gamma - 1.)) 
-u = np.full((N, ), u)
-ds = grp.create_dataset('InternalEnergy', (N,), 'f')
+u = v_c ** 2 / (2.0 * (gamma - 1.0))
+u = np.full((N,), u)
+ds = grp.create_dataset("InternalEnergy", (N,), "f")
 ds[()] = u
 u = np.zeros(1)
 
 # Particle IDs
 ids = 1 + np.linspace(0, N, N, endpoint=False)
-ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
+ds = grp.create_dataset("ParticleIDs", (N,), "L")
 ds[()] = ids
 
 file.close()

examples/Cooling/CoolingHalo/makeIC_random_box.py

 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2016 Stefan Arridge (stefan.arridge@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/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Stefan Arridge (stefan.arridge@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/>.
+#
+##############################################################################
 
 import h5py
 import sys
@@ -24,76 +24,84 @@ import math
 import random
 
 # Generates N particles in a spherically symmetric distribution with density profile ~r^(-2)
-# usage: python makeIC.py 1000: generate 1000 particles
+# usage: python3 makeIC.py 1000: generate 1000 particles
 
 # Some constants
 
-OMEGA = 0.3 # Cosmological matter fraction at z = 0
+OMEGA = 0.3  # Cosmological matter fraction at z = 0
 PARSEC_IN_CGS = 3.0856776e18
 KM_PER_SEC_IN_CGS = 1.0e5
 CONST_G_CGS = 6.672e-8
-h = 0.67777 # hubble parameter
-gamma = 5./3.
+h = 0.67777  # hubble parameter
+gamma = 5.0 / 3.0
 eta = 1.2349
 
 # First set unit velocity and then the circular velocity parameter for the isothermal potential
-const_unit_velocity_in_cgs = 1.e5 #kms^-1
+const_unit_velocity_in_cgs = 1.0e5  # kms^-1
 
-v_c = 200.
+v_c = 200.0
 v_c_cgs = v_c * const_unit_velocity_in_cgs
 
 # Now we use this to get the virial mass and virial radius, which we will set to be the unit mass and radius
 
 # Find H_0, the inverse Hubble time, in cgs
 
-H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+H_0_cgs = 100.0 * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
 
 # From this we can find the virial radius, the radius within which the average density of the halo is
 # 200. * the mean matter density
 
-r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
+r_vir_cgs = v_c_cgs / (10.0 * H_0_cgs * np.sqrt(OMEGA))
 
 # Now get the virial mass
 
-M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+M_vir_cgs = r_vir_cgs * v_c_cgs ** 2 / CONST_G_CGS
 
 # Now set the unit length and mass
 
 const_unit_mass_in_cgs = M_vir_cgs
 const_unit_length_in_cgs = r_vir_cgs
 
-print "UnitMass_in_cgs:     ", const_unit_mass_in_cgs 
-print "UnitLength_in_cgs:   ", const_unit_length_in_cgs
-print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
-
-#derived quantities
-const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
-print "UnitTime_in_cgs:     ", const_unit_time_in_cgs
-const_G                = ((CONST_G_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
-print 'G=', const_G
+print("UnitMass_in_cgs:     ", const_unit_mass_in_cgs)
+print("UnitLength_in_cgs:   ", const_unit_length_in_cgs)
+print("UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs)
+
+# derived quantities
+const_unit_time_in_cgs = const_unit_length_in_cgs / const_unit_velocity_in_cgs
+print("UnitTime_in_cgs:     ", const_unit_time_in_cgs)
+const_G = (
+    CONST_G_CGS
+    * const_unit_mass_in_cgs
+    * const_unit_time_in_cgs
+    * const_unit_time_in_cgs
+    / (const_unit_length_in_cgs * const_unit_length_in_cgs * const_unit_length_in_cgs)
+)
+print("G=", const_G)
 
 # Parameters
-periodic= 1            # 1 For periodic box
-boxSize = 4.          
-G       = const_G 
-N       = int(sys.argv[1])  # Number of particles
+periodic = 1  # 1 For periodic box
+boxSize = 4.0
+G = const_G
+N = int(sys.argv[1])  # Number of particles
 
 # Create the file
 filename = "random_box.hdf5"
-file = h5py.File(filename, 'w')
+file = h5py.File(filename, "w")
 
-#Units
+# Units
 grp = file.create_group("/Units")
 grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
-grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs 
-grp.attrs["Unit time in cgs (U_t)"] = const_unit_length_in_cgs / const_unit_velocity_in_cgs
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs
+grp.attrs["Unit time in cgs (U_t)"] = (
+    const_unit_length_in_cgs / const_unit_velocity_in_cgs
+)
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
 
 # Header
 grp = file.create_group("/Header")
 grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] =  [N ,0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total"] = [N, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_ThisFile"] = [N, 0, 0, 0, 0, 0]
 grp.attrs["Time"] = 0.0
@@ -111,54 +119,54 @@ grp = file.create_group("/PartType0")
 # Positions
 
 # distribute particles randomly in the box
-coords      = np.zeros((N, 3))
-coords[:,0] = boxSize * np.random.rand(N)
-coords[:,1] = boxSize * np.random.rand(N)
-coords[:,2] = boxSize * np.random.rand(N)
-#shift to centre of box
-#coords += np.full((N,3),boxSize/2.)
-print "x range = (%f,%f)" %(np.min(coords[:,0]),np.max(coords[:,0]))
-print "y range = (%f,%f)" %(np.min(coords[:,1]),np.max(coords[:,1]))
-print "z range = (%f,%f)" %(np.min(coords[:,2]),np.max(coords[:,2]))
-
-print np.mean(coords[:,0])
-print np.mean(coords[:,1])
-print np.mean(coords[:,2])
-
-ds = grp.create_dataset('Coordinates', (N, 3), 'd')
+coords = np.zeros((N, 3))
+coords[:, 0] = boxSize * np.random.rand(N)
+coords[:, 1] = boxSize * np.random.rand(N)
+coords[:, 2] = boxSize * np.random.rand(N)
+# shift to centre of box
+# coords += np.full((N,3),boxSize/2.)
+print("x range = (%f,%f)" % (np.min(coords[:, 0]), np.max(coords[:, 0])))
+print("y range = (%f,%f)" % (np.min(coords[:, 1]), np.max(coords[:, 1])))
+print("z range = (%f,%f)" % (np.min(coords[:, 2]), np.max(coords[:, 2])))
+
+print(np.mean(coords[:, 0]))
+print(np.mean(coords[:, 1]))
+print(np.mean(coords[:, 2]))
+
+ds = grp.create_dataset("Coordinates", (N, 3), "d")
 ds[()] = coords
 coords = np.zeros(1)
 
 # All velocities set to zero
-v = np.zeros((N,3))
-ds = grp.create_dataset('Velocities', (N, 3), 'f')
+v = np.zeros((N, 3))
+ds = grp.create_dataset("Velocities", (N, 3), "f")
 ds[()] = v
 v = np.zeros(1)
 
 # All particles of equal mass
-mass = 1. / N
-m = np.full((N,),mass)
-ds = grp.create_dataset('Masses', (N, ), 'f')
+mass = 1.0 / N
+m = np.full((N,), mass)
+ds = grp.create_dataset("Masses", (N,), "f")
 ds[()] = m
 m = np.zeros(1)
 
 # Smoothing lengths
-l = (boxSize**3 / N)**(1./3.) #local mean inter-particle separation
-h = np.full((N, ), eta * l)
-ds = grp.create_dataset('SmoothingLength', (N,), 'f')
+l = (boxSize ** 3 / N) ** (1.0 / 3.0)  # local mean inter-particle separation
+h = np.full((N,), eta * l)
+ds = grp.create_dataset("SmoothingLength", (N,), "f")
 ds[()] = h
 h = np.zeros(1)
 
 # Internal energies
-u = v_c**2 / (2. * (gamma - 1.)) 
-u = np.full((N, ), u)
-ds = grp.create_dataset('InternalEnergy', (N,), 'f')
+u = v_c ** 2 / (2.0 * (gamma - 1.0))
+u = np.full((N,), u)
+ds = grp.create_dataset("InternalEnergy", (N,), "f")
 ds[()] = u
 u = np.zeros(1)
 
 # Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False, dtype='L')
-ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
+ids = 1 + np.linspace(0, N, N, endpoint=False, dtype="L")
+ds = grp.create_dataset("ParticleIDs", (N,), "L")
 ds[()] = ids
 
 file.close()

examples/Cooling/CoolingHalo/run.sh

@@ -2,12 +2,12 @@
 
 # Generate the initial conditions if they are not present.
 echo "Generating initial conditions for the isothermal potential box example..."
-python makeIC.py 10000 
+python3 makeIC.py 10000 
 
-../../swift --external-gravity --hydro --cooling --threads=16 cooling_halo.yml 2>&1 | tee output.log
+../../../swift --external-gravity --hydro --cooling --threads=16 cooling_halo.yml 2>&1 | tee output.log
 
-python radial_profile.py 2. 200 100
+python3 radial_profile.py 2. 200 100
 
-python internal_energy_profile.py 2. 200 100
+python3 internal_energy_profile.py 2. 200 100
 
-python test_energy_conservation.py 2. 200 100
+python3 test_energy_conservation.py 2. 200 100