examples/HydroTests/BlobTest_Braspenning2022/plot_cloud_evolution.py

+from swiftsimio import load
+from swiftsimio.visualisation.volume_render import render_gas
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+import os
+from unyt import s, pc, Myr, kb, mp, cm, K, erg, g, G, km
+import h5py
+import glob
+
+
+def analyse_snapshot(i, chi, rho_cl0, m_cl0, resolution):
+    print("analysing snapshot", i)
+    try:
+        data = load("./{}/blob_{:0=4d}.hdf5".format(resolution, i))
+    except OSError:
+        print("snapshot {} is not found".format(i))
+        exit()
+
+    cloud = data.gas.densities.value > rho_cl0 / 3
+    m_cl = np.sum(data.gas.masses.value[cloud])
+    dense_mass = m_cl / m_cl0
+
+    # Resample onto a 16^3 grid
+    # Naive histogram
+    x_cl = np.median(data.gas.coordinates[:, 0].value[cloud])
+    bins_y = np.linspace(0, 1, 17)
+    bins_z = np.linspace(0, 1, 17)
+    bins_x = np.linspace(0, 4, 4 * 16 + 1)
+    # Shift x_bins such that x_cl is in the centre of a bin
+    bin_space_x = 1 / 16
+    bins_x += x_cl - bin_space_x / 2
+    bins_x = np.sort(bins_x % 4)
+    print(bins_x)
+    binned_masses, _ = np.histogramdd(
+        data.gas.coordinates.v, weights=data.gas.masses.v, bins=[bins_x, bins_y, bins_z]
+    )
+    binned_masses = binned_masses * data.gas.masses.units
+    densities_resampled = binned_masses / (data.metadata.boxsize[1] / 16) ** 3
+    densities_resampled = densities_resampled.to(data.gas.densities.units)
+
+    cloud_resampled = densities_resampled.value > rho_cl0 / 3
+    m_cl_resampled = np.sum(binned_masses.value[cloud_resampled])
+    # For volume rendering we compute the mass from the density and the volume
+    dense_mass_resampled_16 = m_cl_resampled / m_cl0
+
+    # Resample onto a 32^3 grid
+    # Naive histogram
+    x_cl = np.median(data.gas.coordinates[:, 0].value[cloud])
+    bins_y = np.linspace(0, 1, 33)
+    bins_z = np.linspace(0, 1, 33)
+    bins_x = np.linspace(0, 4, 4 * 32 + 1)
+    # Shift x_bins such that x_cl is in the centre of a bin
+    bin_space_x = 1 / 32
+    bins_x += x_cl - bin_space_x / 2
+    bins_x = np.sort(bins_x % 4)
+    print(bins_x)
+    binned_masses, _ = np.histogramdd(
+        data.gas.coordinates.v, weights=data.gas.masses.v, bins=[bins_x, bins_y, bins_z]
+    )
+    binned_masses = binned_masses * data.gas.masses.units
+    densities_resampled = binned_masses / (data.metadata.boxsize[1] / 32) ** 3
+    densities_resampled = densities_resampled.to(data.gas.densities.units)
+
+    cloud_resampled = densities_resampled.value > rho_cl0 / 3
+    m_cl_resampled = np.sum(binned_masses.value[cloud_resampled])
+    # For volume rendering we compute the mass from the density and the volume
+    dense_mass_resampled_32 = m_cl_resampled / m_cl0
+
+    t = data.metadata.time.to(s).value
+
+    # Convert internal energies to temperatures
+    temperatures = (((5 / 3 - 1) * data.gas.internal_energies / kb) * mp).value
+    data.gas.temperatures = temperatures
+
+    if chi == 100:
+        m_t_intermediate = (
+            np.sum(
+                data.gas.masses.value[
+                    (data.gas.temperatures > 10 ** 4.5)
+                    & (data.gas.temperatures < 10 ** 5.5)
+                ]
+            )
+            / m_cl0
+        )
+    elif chi == 10:
+        m_t_intermediate = (
+            np.sum(
+                data.gas.masses.value[
+                    (data.gas.temperatures > 10 ** 4.25)
+                    & (data.gas.temperatures < 10 ** 4.75)
+                ]
+            )
+            / m_cl0
+        )
+
+    """
+    resample the intermediate temperature gas to 32^3 resolution
+    """
+    x_cl = np.median(data.gas.coordinates[:, 0].value[cloud])
+    bins_y = np.linspace(0, 1, 33)
+    bins_z = np.linspace(0, 1, 33)
+    bins_x = np.linspace(0, 4, 4 * 32 + 1)
+    # Shift x_bins such that x_cl is in the centre of a bin
+    bin_space_x = 1 / 32
+    bins_x += x_cl - bin_space_x / 2
+    bins_x = np.sort(bins_x % 4)
+    print(bins_x)
+    binned_masses, _ = np.histogramdd(
+        data.gas.coordinates.v, weights=data.gas.masses.v, bins=[bins_x, bins_y, bins_z]
+    )
+    binned_masstemp, _ = np.histogramdd(
+        data.gas.coordinates.v,
+        weights=data.gas.masses.v * data.gas.temperatures,
+        bins=[bins_x, bins_y, bins_z],
+    )
+    binned_temperatures = binned_masstemp / binned_masses
+    binned_masses = binned_masses * data.gas.masses.units
+
+    if chi == 100:
+        mask = (binned_temperatures > 10 ** 4.5) & (binned_temperatures < 10 ** 5.5)
+        m_t_intermediate_resampled_32 = np.sum(binned_masses.v[mask]) / m_cl0
+    elif chi == 10:
+        mask = (binned_temperatures > 10 ** 4.25) & (binned_temperatures < 10 ** 4.75)
+        m_t_intermediate_resampled_32 = np.sum(binned_masses.v[mask]) / m_cl0
+
+    """
+    resample the intermediate temperature gas to 16^3 resolution 
+    """
+    x_cl = np.median(data.gas.coordinates[:, 0].value[cloud])
+    bins_y = np.linspace(0, 1, 17)
+    bins_z = np.linspace(0, 1, 17)
+    bins_x = np.linspace(0, 4, 4 * 16 + 1)
+    # Shift x_bins such that x_cl is in the centre of a bin
+    bin_space_x = 1 / 16
+    bins_x += x_cl - bin_space_x / 2
+    bins_x = np.sort(bins_x % 4)
+    print(bins_x)
+    binned_masses, _ = np.histogramdd(
+        data.gas.coordinates.v, weights=data.gas.masses.v, bins=[bins_x, bins_y, bins_z]
+    )
+    binned_masstemp, _ = np.histogramdd(
+        data.gas.coordinates.v,
+        weights=data.gas.masses.v * data.gas.temperatures,
+        bins=[bins_x, bins_y, bins_z],
+    )
+    binned_temperatures = binned_masstemp / binned_masses
+    binned_masses = binned_masses * data.gas.masses.units
+
+    if chi == 100:
+        mask = (binned_temperatures > 10 ** 4.5) & (binned_temperatures < 10 ** 5.5)
+        m_t_intermediate_resampled_16 = np.sum(binned_masses.v[mask]) / m_cl0
+    elif chi == 10:
+        mask = (binned_temperatures > 10 ** 4.25) & (binned_temperatures < 10 ** 4.75)
+        m_t_intermediate_resampled_16 = np.sum(binned_masses.v[mask]) / m_cl0
+
+    return (
+        t,
+        dense_mass,
+        m_t_intermediate,
+        dense_mass_resampled_16,
+        m_t_intermediate_resampled_16,
+        dense_mass_resampled_32,
+        m_t_intermediate_resampled_32,
+    )
+
+
+def load_first_snapshot(chi, resolution):
+    # Load snapshot 0
+    try:
+        data_0 = load("./{}/blob_0000.hdf5".format(resolution))
+    except OSError:
+        print("snapshot 0 is not found")
+
+    v_wind0 = data_0.gas.velocities[
+        0, 0
+    ]  # All wind particles have the same initial velocity
+    r_cl0 = data_0.metadata.boxsize[1] * 0.1  # The cloud is 0.1 * boxsize
+
+    rho_cl0 = chi * np.percentile(data_0.gas.densities.value, 50)
+    cloud_0 = (
+        data_0.gas.densities.value > rho_cl0 / 3
+    )  # Create a mask for the dense gas
+    m_cl0 = np.sum(data_0.gas.masses.value[cloud_0])
+
+    t_cc = (np.sqrt(chi) * r_cl0 / v_wind0).to(s)  # Cloud crushing time
+
+    # This will give the number of snapshots + 1 IC file
+    hdf5counter = len(glob.glob1("./{}".format(resolution), "*.hdf5"))
+    # Subtract the IC file from the count
+    steps = hdf5counter - 1
+
+    return rho_cl0, m_cl0, steps, t_cc
+
+
+# What is the number of (particles)^1/3?
+resolution = 128
+# What is the density contrast?
+chi = 10
+
+rho_cl0, m_cl0, steps, t_cc = load_first_snapshot(chi, resolution)
+
+time = np.zeros(steps)
+dense_mass = np.zeros(steps)
+intermediate_temp_mass = np.zeros(steps)
+
+dense_mass_16 = np.zeros(steps)
+intermediate_temp_mass_16 = np.zeros(steps)
+
+dense_mass_32 = np.zeros(steps)
+intermediate_temp_mass_32 = np.zeros(steps)
+
+for i in range(steps):
+    time[i], dense_mass[i], intermediate_temp_mass[i], dense_mass_16[
+        i
+    ], intermediate_temp_mass_16[i], dense_mass_32[i], intermediate_temp_mass_32[
+        i
+    ] = analyse_snapshot(
+        i, chi, rho_cl0, m_cl0, resolution
+    )
+
+fig, ax = plt.subplots(3, 2, figsize=(10, 15))
+
+ax[0, 0].plot(time / t_cc, dense_mass)
+ax[0, 0].set_title("Dense Mass")
+
+ax[0, 1].plot(time / t_cc, intermediate_temp_mass)
+ax[0, 1].set_title("Intermediate-temperature Mass")
+
+ax[1, 0].plot(time / t_cc, dense_mass_16)
+ax[1, 0].set_title("Dense Mass resampled 16^3")
+
+ax[1, 1].plot(time / t_cc, intermediate_temp_mass_16)
+ax[1, 1].set_title("Intermediate-temperature Mass resampled 16^3")
+
+ax[2, 0].plot(time / t_cc, dense_mass_32)
+ax[2, 0].set_title("Dense Mass resampled 32^3")
+
+ax[2, 1].plot(time / t_cc, intermediate_temp_mass_32)
+ax[2, 1].set_title("Intermediate-temperature Mass resampled 32^3")
+
+plt.savefig("cloud_evolution")
+plt.close()

examples/HydroTests/Diffusion_1D/run.sh

@@ -4,8 +4,8 @@
 if [ ! -e diffusion.hdf5 ]
 then
     echo "Generating initial conditions for the Sedov blast example..."
-    python makeIC.py
+    python3 makeIC.py
 fi
 
 # Run SWIFT
-../../swift --hydro --limiter --threads=1 diffusion.yml 2>&1 | tee output.log
+../../../swift --hydro --limiter --threads=1 diffusion.yml 2>&1 | tee output.log

examples/HydroTests/Diffusion_1D/run_set.sh

 #!/bin/bash
 
 beta=(1.0 0.1 0.01 0.001)
-swift_location="../../swift"
+swift_location="../../../swift"
 flags="--hydro --limiter --threads=2"
 parameter="diffusion.yml"
 parameter_fixed="diffusion_fixed_alpha.yml"
@@ -13,12 +13,12 @@ do
     mkdir beta_$i
     cd beta_$i
 
-    python ../$make_ic_script
+    python3 ../$make_ic_script
 
     ../$swift_location $flags -P "SPH:diffusion_beta:${i}" ../$parameter
-    python ../$plot_script
+    python3 ../$plot_script
     ../$swift_location $flags -P "SPH:diffusion_beta:${i}" ../$parameter_fixed
-    python ../$plot_script -s diffusion_fixed_alpha -o diffusion_fixed_alpha.png
+    python3 ../$plot_script -s diffusion_fixed_alpha -o diffusion_fixed_alpha.png
 
     rm *.hdf5
 

examples/HydroTests/EvrardCollapse_3D/getReference.sh

 #! /bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/evrardCollapse3D_exact.txt
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/evrardCollapse3D_exact.txt

examples/HydroTests/EvrardCollapse_3D/makeIC.py

@@ -32,36 +32,36 @@ parser.add_argument(
          Number of particles to be used in the Evrard collapse.
          """,
     required=False,
-    default=100000
+    default=100000,
 )
 
 args = vars(parser.parse_args())
 
 # Parameters
-gamma = 5. / 3.      # Gas adiabatic index
-M = 1.  # total mass of the sphere
-R = 1.               # radius of the sphere
-u0 = 0.05 / M        # initial thermal energy
-fileName = "evrard.hdf5" 
+gamma = 5.0 / 3.0  # Gas adiabatic index
+M = 1.0  # total mass of the sphere
+R = 1.0  # radius of the sphere
+u0 = 0.05 / M  # initial thermal energy
+fileName = "evrard.hdf5"
 numPart = int(args["nparts"])
 
 r = R * sqrt(random.random(numPart))
-phi = 2. * pi * random.random(numPart)
-cos_theta = 2. * random.random(numPart) - 1.
+phi = 2.0 * pi * random.random(numPart)
+cos_theta = 2.0 * random.random(numPart) - 1.0
 
-sin_theta = sqrt(1. - cos_theta**2)
+sin_theta = sqrt(1.0 - cos_theta ** 2)
 cos_phi = cos(phi)
 sin_phi = sin(phi)
 
 pos = zeros((numPart, 3))
-pos[:,0] = r * sin_theta * cos_phi
-pos[:,1] = r * sin_theta * sin_phi
-pos[:,2] = r * cos_theta
+pos[:, 0] = r * sin_theta * cos_phi
+pos[:, 1] = r * sin_theta * sin_phi
+pos[:, 2] = r * cos_theta
 
 # shift particles to put the sphere in the centre of the box
-pos += array([50. * R, 50. * R, 50. * R])
+pos += array([50.0 * R, 50.0 * R, 50.0 * R])
 
-h = ones(numPart) * 2. * R / numPart**(1. / 3.)
+h = ones(numPart) * 2.0 * R / numPart ** (1.0 / 3.0)
 
 # Generate extra arrays
 v = zeros((numPart, 3))
@@ -69,15 +69,15 @@ ids = linspace(1, numPart, numPart)
 m = ones(numPart) * M / numPart
 u = ones(numPart) * u0
 
-#--------------------------------------------------
+# --------------------------------------------------
 
-#File
-file = h5py.File(fileName, 'w')
+# File
+file = h5py.File(fileName, "w")
 
 # Header
 grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = [100. * R, 100. * R, 100. * R]
-grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["BoxSize"] = [100.0 * R, 100.0 * R, 100.0 * R]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
 grp.attrs["Time"] = 0.0
@@ -86,21 +86,21 @@ grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 grp.attrs["Flag_Entropy_ICs"] = 0
 grp.attrs["Dimension"] = 3
 
-#Units
+# Units
 grp = file.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 1.
-grp.attrs["Unit mass in cgs (U_M)"] = 1.
-grp.attrs["Unit time in cgs (U_t)"] = 1.
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+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
+# Particle group
 grp = file.create_group("/PartType0")
-grp.create_dataset('Coordinates', data=pos, dtype='d')
-grp.create_dataset('Velocities', data=v, dtype='f')
-grp.create_dataset('Masses', data=m, dtype='f')
-grp.create_dataset('SmoothingLength', data=h, dtype='f')
-grp.create_dataset('InternalEnergy', data=u, dtype='f')
-grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+grp.create_dataset("Coordinates", data=pos, dtype="d")
+grp.create_dataset("Velocities", data=v, dtype="f")
+grp.create_dataset("Masses", data=m, dtype="f")
+grp.create_dataset("SmoothingLength", data=h, dtype="f")
+grp.create_dataset("InternalEnergy", data=u, dtype="f")
+grp.create_dataset("ParticleIDs", data=ids, dtype="L")
 
 file.close()

examples/Logger/SimpleOrbits/makeIC.py → examples/HydroTests/EvrardCollapse_3D/plotEnergy.py

 ###############################################################################
 # This file is part of SWIFT.
-# Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
+# Copyright (c) 2023 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
@@ -15,61 +15,44 @@
 # 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 numpy as np
-from astropy import units
-from swiftsimio import Writer
-import unyt
-
-np.random.seed(50)
+import matplotlib
 
-# parameters
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy import stats
+import sys
 
-filename = "simple_orbits.hdf5"
-num_part = 1
-masses = 1.
-# If changed, need to update simple_orbits.yml
-M = units.solMass.to("earthMass")
-M = float(M)
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
-min_r = 0.2
-max_r = 5
-boxsize = 2 * max_r
+data = np.loadtxt("statistics.txt")
+t = data[:, 1]
+E_kin = data[:, 13]
+E_int = data[:, 14]
+E_pot = data[:, 15]
 
-u_l = 1.49e13  # AU
-u_m = 5.97e27  # Earth mass
-u_v = 474047.  # AU / yr
-u_t = u_l / u_v
-G = 1.189972e-04  # grav. const.
+E_tot = E_kin + E_int + E_pot
 
-# generate the coordinates
-dist = np.random.rand(num_part) * (max_r - min_r) + min_r
-angle = np.random.rand(num_part) * 2 * np.pi
-# more easy to do it in 2D, therefore coords[:, 2] == 0
-coords = np.zeros((num_part, 3))
-coords[:, 0] = dist * np.sin(angle)
-coords[:, 1] = dist * np.cos(angle)
-coords += boxsize * 0.5
+plt.figure(figsize=(5, 5))
 
-# generate the masses
-m = np.ones(num_part) * masses
+plt.subplot(211)
+plt.plot(t, E_tot, "k", label="total")
+plt.plot(t, E_int, label="internal")
+plt.plot(t, E_kin, label="kinetic")
+plt.plot(t, E_pot, label="potential")
 
-# generate the velocities
-sign = np.random.rand(num_part)
-sign[sign < 0.5] = -1
-sign[sign >= 0.5] = 1
+plt.legend(loc="upper left")
 
-v = np.zeros((num_part, 3))
-v[:, 0] = sign * np.sqrt(G * M / (dist * (1 + np.tan(angle)**2)))
-v[:, 1] = - np.tan(angle) * v[:, 0]
+plt.ylabel("Energy")
+plt.xlabel("Time")
 
-# Write the snapshot
-units = unyt.UnitSystem("Planets", unyt.AU, unyt.mearth, unyt.yr)
+plt.subplot(212)
+plt.plot(t, E_tot / E_tot[0], "k")
 
-snapshot = Writer(units, boxsize * unyt.AU)
-snapshot.dark_matter.coordinates = coords * unyt.AU
-snapshot.dark_matter.velocities = v * unyt.AU / unyt.yr
-snapshot.dark_matter.masses = m * unyt.mearth
+plt.ylim(0.98, 1.02)
+plt.ylabel("Total energy w.r.t. t=0")
+plt.xlabel("Time")
 
-snapshot.write(filename)
+plt.savefig("Energy.png")

examples/HydroTests/EvrardCollapse_3D/plotSolution.py

 ###############################################################################
 # This file is part of SWIFT.
-# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+# Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
 #               2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
-# 
+#
 # 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/>.
-# 
+#
 ################################################################################
 
 # Compares the swift result for the 2D spherical Sod shock with a high
 # resolution 2D reference result
 
 import matplotlib
+
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 from scipy import stats
 import h5py
+import sys
 
 # Parameters
-gas_gamma = 5./3.      # Polytropic index
-rho_L = 1.             # Density left state
-rho_R = 0.125          # Density right state
-v_L = 0.               # Velocity left state
-v_R = 0.               # Velocity right state
-P_L = 1.               # Pressure left state
-P_R = 0.1              # Pressure right state
-
-# 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' : (9.90,6.45),
-'figure.subplot.left'    : 0.045,
-'figure.subplot.right'   : 0.99,
-'figure.subplot.bottom'  : 0.05,
-'figure.subplot.top'     : 0.99,
-'figure.subplot.wspace'  : 0.15,
-'figure.subplot.hspace'  : 0.12,
-'lines.markersize' : 6,
-'lines.linewidth' : 3.,
-'text.latex.unicode': True
-}
-rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+gas_gamma = 5.0 / 3.0  # Polytropic index
+rho_L = 1.0  # Density left state
+rho_R = 0.125  # Density right state
+v_L = 0.0  # Velocity left state
+v_R = 0.0  # Velocity right state
+P_L = 1.0  # Pressure left state
+P_R = 0.1  # Pressure right state
+
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
 # Read the simulation data
-sim = h5py.File("evrard_%04d.hdf5"%snap, "r")
+sim = h5py.File("evrard_%04d.hdf5" % snap, "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
 time = sim["/Header"].attrs["Time"][0]
 scheme = sim["/HydroScheme"].attrs["Scheme"]
@@ -71,102 +54,148 @@ eta = sim["/HydroScheme"].attrs["Kernel eta"]
 git = sim["Code"].attrs["Git Revision"]
 
 coords = sim["/PartType0/Coordinates"]
-x = sqrt((coords[:,0] - 0.5 * boxSize)**2 + (coords[:,1] - 0.5 * boxSize)**2 + \
-         (coords[:,2] - 0.5 * boxSize)**2)
+x = np.sqrt(
+    (coords[:, 0] - 0.5 * boxSize) ** 2
+    + (coords[:, 1] - 0.5 * boxSize) ** 2
+    + (coords[:, 2] - 0.5 * boxSize) ** 2
+)
 vels = sim["/PartType0/Velocities"]
-v = sqrt(vels[:,0]**2 + vels[:,1]**2 + vels[:,2]**2)
+v = np.sqrt(vels[:, 0] ** 2 + vels[:, 1] ** 2 + vels[:, 2] ** 2)
 u = sim["/PartType0/InternalEnergies"][:]
 S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
 rho = sim["/PartType0/Densities"][:]
 
 # Bin the data
-x_bin_edge = logspace(-3., log10(2.), 100)
-x_bin = 0.5*(x_bin_edge[1:] + x_bin_edge[:-1])
-rho_bin,_,_ = stats.binned_statistic(x, rho, statistic='mean', bins=x_bin_edge)
-v_bin,_,_ = stats.binned_statistic(x, v, statistic='mean', bins=x_bin_edge)
-P_bin,_,_ = stats.binned_statistic(x, P, statistic='mean', bins=x_bin_edge)
-S_bin,_,_ = stats.binned_statistic(x, S, statistic='mean', bins=x_bin_edge)
-u_bin,_,_ = stats.binned_statistic(x, u, statistic='mean', bins=x_bin_edge)
-rho2_bin,_,_ = stats.binned_statistic(x, rho**2, statistic='mean', bins=x_bin_edge)
-v2_bin,_,_ = stats.binned_statistic(x, v**2, statistic='mean', bins=x_bin_edge)
-P2_bin,_,_ = stats.binned_statistic(x, P**2, statistic='mean', bins=x_bin_edge)
-S2_bin,_,_ = stats.binned_statistic(x, S**2, statistic='mean', bins=x_bin_edge)
-u2_bin,_,_ = stats.binned_statistic(x, u**2, statistic='mean', bins=x_bin_edge)
-rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
-v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
-P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
-S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
-u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
-
-ref = loadtxt("evrardCollapse3D_exact.txt")
+x_bin_edge = np.logspace(-3.0, np.log10(2.0), 100)
+x_bin = 0.5 * (x_bin_edge[1:] + x_bin_edge[:-1])
+rho_bin, _, _ = stats.binned_statistic(x, rho, statistic="mean", bins=x_bin_edge)
+v_bin, _, _ = stats.binned_statistic(x, v, statistic="mean", bins=x_bin_edge)
+P_bin, _, _ = stats.binned_statistic(x, P, statistic="mean", bins=x_bin_edge)
+S_bin, _, _ = stats.binned_statistic(x, S, statistic="mean", bins=x_bin_edge)
+u_bin, _, _ = stats.binned_statistic(x, u, statistic="mean", bins=x_bin_edge)
+rho2_bin, _, _ = stats.binned_statistic(x, rho ** 2, statistic="mean", bins=x_bin_edge)
+v2_bin, _, _ = stats.binned_statistic(x, v ** 2, statistic="mean", bins=x_bin_edge)
+P2_bin, _, _ = stats.binned_statistic(x, P ** 2, statistic="mean", bins=x_bin_edge)
+S2_bin, _, _ = stats.binned_statistic(x, S ** 2, statistic="mean", bins=x_bin_edge)
+u2_bin, _, _ = stats.binned_statistic(x, u ** 2, statistic="mean", bins=x_bin_edge)
+rho_sigma_bin = np.sqrt(rho2_bin - rho_bin ** 2)
+v_sigma_bin = np.sqrt(v2_bin - v_bin ** 2)
+P_sigma_bin = np.sqrt(P2_bin - P_bin ** 2)
+S_sigma_bin = np.sqrt(S2_bin - S_bin ** 2)
+u_sigma_bin = np.sqrt(u2_bin - u_bin ** 2)
+
+ref = np.loadtxt("evrardCollapse3D_exact.txt")
 
 # Plot the interesting quantities
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=1,
+    markeredgecolor="none",
+    alpha=0.2,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
+
+errorbar_props = dict(color=binned_color, ms=binned_marker_size, fmt=".", lw=1.2)
 
 # Velocity profile --------------------------------
-subplot(231)
-semilogx(x, -v, '.', color='r', ms=0.2)
-semilogx(ref[:,0], ref[:,2], "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, -v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
-xlim(1.e-3, 2.)
-ylim(-1.7, 0.1)
+plt.subplot(231)
+plt.semilogx(x, -v, **scatter_props)
+plt.semilogx(ref[:, 0], ref[:, 2], "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, -v_bin, yerr=v_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
+plt.xlim(1.0e-3, 2.0)
+plt.ylim(-1.7, 0.1)
 
 # Density profile --------------------------------
-subplot(232)
-loglog(x, rho, '.', color='r', ms=0.2)
-loglog(ref[:,0], ref[:,1], "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
-xlim(1.e-3, 2.)
-ylim(1.e-2, 1.e4)
+plt.subplot(232)
+plt.loglog(x, rho, **scatter_props)
+plt.loglog(ref[:, 0], ref[:, 1], "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.xlim(1.0e-3, 2.0)
+plt.ylim(1.0e-2, 1.0e4)
 
 # Pressure profile --------------------------------
-subplot(233)
-loglog(x, P, '.', color='r', ms=0.2)
-loglog(ref[:,0], ref[:,3], "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Pressure}}~P$", labelpad=0)
-xlim(1.e-3, 2.)
-ylim(1.e-4, 1.e3)
+plt.subplot(233)
+plt.loglog(x, P, **scatter_props)
+plt.loglog(ref[:, 0], ref[:, 3], "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(x_bin, P_bin, yerr=P_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.xlim(1.0e-3, 2.0)
+plt.ylim(1.0e-4, 1.0e3)
 
 # Internal energy profile -------------------------
-subplot(234)
-loglog(x, u, '.', color='r', ms=0.2)
-loglog(ref[:,0], ref[:,3] / ref[:,1] / (gas_gamma - 1.), "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
-xlim(1.e-3, 2.)
-ylim(1.e-2, 2.)
+plt.subplot(234)
+plt.loglog(x, u, **scatter_props)
+plt.loglog(
+    ref[:, 0],
+    ref[:, 3] / ref[:, 1] / (gas_gamma - 1.0),
+    "--",
+    color=line_color,
+    alpha=0.8,
+    lw=1.2,
+)
+plt.errorbar(x_bin, u_bin, yerr=u_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.xlim(1.0e-3, 2.0)
+plt.ylim(1.0e-2, 2.0)
 
 # Entropy profile ---------------------------------
-subplot(235)
-semilogx(x, S, '.', color='r', ms=0.2)
-semilogx(ref[:,0], ref[:,3] / ref[:,1]**gas_gamma, "k--", alpha=0.8, lw=1.2)
-errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=0)
-xlim(1.e-3, 2.)
-ylim(0., 0.25)
+plt.subplot(235)
+plt.semilogx(x, S, **scatter_props)
+plt.semilogx(
+    ref[:, 0],
+    ref[:, 3] / ref[:, 1] ** gas_gamma,
+    "--",
+    color=line_color,
+    alpha=0.8,
+    lw=1.2,
+)
+plt.errorbar(x_bin, S_bin, yerr=S_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+plt.xlim(1.0e-3, 2.0)
+plt.ylim(0.0, 0.25)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
-
-text(-0.49, 0.9, "Evrard collapse with $\\gamma=%.3f$ in 3D\nat $t=%.2f$"%(gas_gamma,time), fontsize=10)
-plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-tight_layout()
-savefig("EvrardCollapse.png", dpi=200)
+plt.subplot(236, frameon=False)
+
+text_fontsize = 5
+
+plt.text(
+    -0.45,
+    0.9,
+    "Evrard collapse with $\\gamma=%.3f$ in 3D\nat $t=%.2f$" % (gas_gamma, time),
+    fontsize=text_fontsize,
+)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
+)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("EvrardCollapse.png", dpi=200)

examples/HydroTests/EvrardCollapse_3D/run.sh

@@ -4,11 +4,11 @@
 if [ ! -e evrard.hdf5 ]
 then
     echo "Generating initial conditions for the Evrard collapse example..."
-    python makeIC.py
+    python3 makeIC.py
 fi
 
 # Run SWIFT
-../../swift --hydro --self-gravity --threads=4 evrard.yml 2>&1 | tee output.log
+../../../swift --hydro --self-gravity --threads=4 evrard.yml 2>&1 | tee output.log
 
 # Get the high resolution 1D reference result if not present.
 if [ ! -e evrardCollapse3D_exact.txt ]
@@ -18,4 +18,4 @@ then
 fi
 
 # Plot the solution
-python plotSolution.py 8
+python3 plotSolution.py 8

examples/HydroTests/Gradients/makeICs.py

@@ -26,7 +26,7 @@ import sys
 # reconstruction
 
 # Parameters
-gamma = 5./3.     # Gas adiabatic index
+gamma = 5.0 / 3.0  # Gas adiabatic index
 gridtype = "cartesian"
 if len(sys.argv) > 1:
     gridtype = sys.argv[1]
@@ -35,116 +35,116 @@ if len(sys.argv) > 1:
 if gridtype == "stretched":
     fileName = "Gradients_stretched.hdf5"
     factor = 8
-    boxSize = [ 1.0 , 1.0/factor , 1.0/factor ]
+    boxSize = [1.0, 1.0 / factor, 1.0 / factor]
     L = 20
-    nx1 = factor*L/2
+    nx1 = factor * L / 2
     ny1 = L
     nz1 = L
-    numfac = 2.
-    nx2 = int(nx1/numfac)
-    ny2 = int(ny1/numfac)
-    nz2 = int(nz1/numfac)
-    npart = nx1*ny1*nz1 + nx2*ny2*nz2
+    numfac = 2.0
+    nx2 = int(nx1 / numfac)
+    ny2 = int(ny1 / numfac)
+    nz2 = int(nz1 / numfac)
+    npart = nx1 * ny1 * nz1 + nx2 * ny2 * nz2
     vol = boxSize[0] * boxSize[1] * boxSize[2]
-    partVol1 = 0.5*vol/(nx1*ny1*nz1)
-    partVol2 = 0.5*vol/(nx2*ny2*nz2)
+    partVol1 = 0.5 * vol / (nx1 * ny1 * nz1)
+    partVol2 = 0.5 * vol / (nx2 * ny2 * nz2)
 
-    coords = np.zeros((npart,3))
-    h = np.zeros((npart,1))
-    ids = np.zeros((npart,1), dtype='L')
+    coords = np.zeros((npart, 3))
+    h = np.zeros((npart, 1))
+    ids = np.zeros((npart, 1), dtype="L")
     idx = 0
-    dcell = 0.5/nx1
+    dcell = 0.5 / nx1
     for i in range(nx1):
         for j in range(ny1):
             for k in range(nz1):
-                coords[idx,0] = (i+0.5)*dcell
-                coords[idx,1] = (j+0.5)*dcell
-                coords[idx,2] = (k+0.5)*dcell
-                h[idx] = 0.56/nx1
+                coords[idx, 0] = (i + 0.5) * dcell
+                coords[idx, 1] = (j + 0.5) * dcell
+                coords[idx, 2] = (k + 0.5) * dcell
+                h[idx] = 0.56 / nx1
                 ids[idx] = idx
                 idx += 1
-    dcell = 0.5/nx2
+    dcell = 0.5 / nx2
     for i in range(nx2):
         for j in range(ny2):
             for k in range(nz2):
-                coords[idx,0] = 0.5+(i+0.5)*dcell
-                coords[idx,1] = (j+0.5)*dcell
-                coords[idx,2] = (k+0.5)*dcell
-                h[idx] = 0.56/nx2
+                coords[idx, 0] = 0.5 + (i + 0.5) * dcell
+                coords[idx, 1] = (j + 0.5) * dcell
+                coords[idx, 2] = (k + 0.5) * dcell
+                h[idx] = 0.56 / nx2
                 ids[idx] = idx
                 idx += 1
 
 # cartesian box ################################################################
 if gridtype == "cartesian":
     fileName = "Gradients_cartesian.hdf5"
-    boxSize = [ 1.0 , 1.0 , 1.0 ]
+    boxSize = [1.0, 1.0, 1.0]
     nx = 20
-    npart = nx**3
-    partVol = 1./npart
-    coords = np.zeros((npart,3))
-    h = np.zeros((npart,1))
-    ids = np.zeros((npart,1), dtype='L')
+    npart = nx ** 3
+    partVol = 1.0 / npart
+    coords = np.zeros((npart, 3))
+    h = np.zeros((npart, 1))
+    ids = np.zeros((npart, 1), dtype="L")
     idx = 0
-    dcell = 1./nx
+    dcell = 1.0 / nx
     for i in range(nx):
         for j in range(nx):
             for k in range(nx):
-                coords[idx,0] = (i+0.5)*dcell
-                coords[idx,1] = (j+0.5)*dcell
-                coords[idx,2] = (k+0.5)*dcell
-                h[idx] = 1.12/nx
+                coords[idx, 0] = (i + 0.5) * dcell
+                coords[idx, 1] = (j + 0.5) * dcell
+                coords[idx, 2] = (k + 0.5) * dcell
+                h[idx] = 1.12 / nx
                 ids[idx] = idx
                 idx += 1
 
 # random box ###################################################################
 if gridtype == "random":
     fileName = "Gradients_random.hdf5"
-    boxSize = [ 1.0 , 1.0 , 1.0 ]
+    boxSize = [1.0, 1.0, 1.0]
     glass = h5py.File("../Glass/glass_50000.hdf5", "r")
     coords = np.array(glass["/PartType0/Coordinates"])
     npart = len(coords)
-    partVol = 1./npart
-    h = np.zeros((npart,1))
-    ids = np.zeros((npart,1), dtype='L')
+    partVol = 1.0 / npart
+    h = np.zeros((npart, 1))
+    ids = np.zeros((npart, 1), dtype="L")
     for i in range(npart):
         h[i] = 0.019
         ids[i] = i
 
-v = np.zeros((npart,3))
-m = np.zeros((npart,1))
-rho = np.zeros((npart,1))
-u = np.zeros((npart,1))
+v = np.zeros((npart, 3))
+m = np.zeros((npart, 1))
+rho = np.zeros((npart, 1))
+u = np.zeros((npart, 1))
 
 for i in range(npart):
-    rhox = coords[i,0]
-    if coords[i,0] < 0.75:
+    rhox = coords[i, 0]
+    if coords[i, 0] < 0.75:
         rhox = 0.75
-    if coords[i,0] < 0.25:
-        rhox = 1.-coords[i,0]
-    rhoy = 1.+boxSize[1]-coords[i,1]
-    if coords[i,1] < 0.75*boxSize[1]:
-        rhoy = 1. + 0.25*boxSize[1]
-    if coords[i,1] < 0.25*boxSize[1]:
-        rhoy = 1.+coords[i,1]
-    rhoz = 1.
+    if coords[i, 0] < 0.25:
+        rhox = 1.0 - coords[i, 0]
+    rhoy = 1.0 + boxSize[1] - coords[i, 1]
+    if coords[i, 1] < 0.75 * boxSize[1]:
+        rhoy = 1.0 + 0.25 * boxSize[1]
+    if coords[i, 1] < 0.25 * boxSize[1]:
+        rhoy = 1.0 + coords[i, 1]
+    rhoz = 1.0
     rho[i] = rhox + rhoy + rhoz
-    P = 1.
-    u[i] = P / ((gamma-1.)*rho[i])
+    P = 1.0
+    u[i] = P / ((gamma - 1.0) * rho[i])
     if gridtype == "stretched":
-        if coords[i,0] < 0.5:
+        if coords[i, 0] < 0.5:
             m[i] = rho[i] * partVol1
         else:
             m[i] = rho[i] * partVol2
     else:
         m[i] = rho[i] * partVol
 
-#File
-file = h5py.File(fileName, 'w')
+# File
+file = h5py.File(fileName, "w")
 
 # Header
 grp = file.create_group("/Header")
 grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] =  [npart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total"] = [npart, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_ThisFile"] = [npart, 0, 0, 0, 0, 0]
 grp.attrs["Time"] = 0.0
@@ -152,21 +152,21 @@ 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, 0, 0, 0, 0, 0]
 
-#Particle group
+# Particle group
 grp = file.create_group("/PartType0")
-ds = grp.create_dataset('Coordinates', (npart, 3), 'd')
+ds = grp.create_dataset("Coordinates", (npart, 3), "d")
 ds[()] = coords
-ds = grp.create_dataset('Velocities', (npart, 3), 'f')
+ds = grp.create_dataset("Velocities", (npart, 3), "f")
 ds[()] = v
-ds = grp.create_dataset('Masses', (npart,1), 'f')
+ds = grp.create_dataset("Masses", (npart, 1), "f")
 ds[()] = m
-ds = grp.create_dataset('Density', (npart,1), 'd')
+ds = grp.create_dataset("Density", (npart, 1), "d")
 ds[()] = rho
-ds = grp.create_dataset('SmoothingLength', (npart,1), 'f')
+ds = grp.create_dataset("SmoothingLength", (npart, 1), "f")
 ds[()] = h
-ds = grp.create_dataset('InternalEnergy', (npart,1), 'd')
+ds = grp.create_dataset("InternalEnergy", (npart, 1), "d")
 ds[()] = u
-ds = grp.create_dataset('ParticleIDs', (npart,1), 'L')
+ds = grp.create_dataset("ParticleIDs", (npart, 1), "L")
 ds[()] = ids[:]
 
 file.close()

examples/HydroTests/Gradients/plot.py

@@ -34,15 +34,15 @@ rho = np.array(f["/PartType0/Densities"])
 gradrho = np.array(f["/PartType0/GradDensity"])
 coords = np.array(f["/PartType0/Coordinates"])
 
-fig, ax = pl.subplots(1,2, sharey=True)
+fig, ax = pl.subplots(1, 2, sharey=True)
 
-ax[0].plot(coords[:,0], rho, "r.", label="density")
-ax[0].plot(coords[:,0], gradrho[:,0], "b.", label="grad density x")
+ax[0].plot(coords[:, 0], rho, "r.", label="density")
+ax[0].plot(coords[:, 0], gradrho[:, 0], "b.", label="grad density x")
 ax[0].set_xlabel("x")
 ax[0].legend(loc="best")
 
-ax[1].plot(coords[:,1], rho, "r.", label="density")
-ax[1].plot(coords[:,1], gradrho[:,1], "b.", label="grad density y")
+ax[1].plot(coords[:, 1], rho, "r.", label="density")
+ax[1].plot(coords[:, 1], gradrho[:, 1], "b.", label="grad density y")
 ax[1].set_xlabel("y")
 ax[1].legend(loc="best")
 

examples/HydroTests/Gradients/run.sh

 #! /bin/bash
 
-python makeICs.py stretched
-../../swift --hydro --threads=2 gradientsStretched.yml
-python plot.py gradients_stretched_0001.hdf5 stretched
+python3 makeICs.py stretched
+../../../swift --hydro --threads=2 gradientsStretched.yml
+python3 plot.py gradients_stretched_0001.hdf5 stretched
 
-python makeICs.py cartesian
-../../swift --hydro --threads=2 gradientsCartesian.yml
-python plot.py gradients_cartesian_0001.hdf5 cartesian
+python3 makeICs.py cartesian
+../../../swift --hydro --threads=2 gradientsCartesian.yml
+python3 plot.py gradients_cartesian_0001.hdf5 cartesian
 
-python makeICs.py random
-../../swift --hydro --threads=2 gradientsRandom.yml
-python plot.py gradients_random_0001.hdf5 random
+python3 makeICs.py random
+../../../swift --hydro --threads=2 gradientsRandom.yml
+python3 plot.py gradients_random_0001.hdf5 random

examples/HydroTests/GreshoVortex_2D/getGlass.sh

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

examples/HydroTests/GreshoVortex_2D/makeIC.py

 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2016 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/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 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 h5py
 from numpy import *
@@ -23,16 +23,16 @@ from numpy import *
 # Generates a swift IC file for the Gresho-Chan vortex in a periodic box
 
 # Parameters
-gamma = 5./3.     # Gas adiabatic index
-rho0 = 1           # Gas density
-P0 = 0.           # Constant additional pressure (should have no impact on the dynamics)
-fileOutputName = "greshoVortex.hdf5" 
+gamma = 5.0 / 3.0  # Gas adiabatic index
+rho0 = 1  # Gas density
+P0 = 0.0  # Constant additional pressure (should have no impact on the dynamics)
+fileOutputName = "greshoVortex.hdf5"
 fileGlass = "glassPlane_128.hdf5"
-#---------------------------------------------------
+# ---------------------------------------------------
 
 # Get position and smoothing lengths from the glass
-fileInput = h5py.File(fileGlass, 'r')
-coords = fileInput["/PartType0/Coordinates"][:,:]
+fileInput = h5py.File(fileGlass, "r")
+coords = fileInput["/PartType0/Coordinates"][:, :]
 h = fileInput["/PartType0/SmoothingLength"][:]
 ids = fileInput["/PartType0/ParticleIDs"][:]
 boxSize = fileInput["/Header"].attrs["BoxSize"][0]
@@ -40,47 +40,46 @@ numPart = size(h)
 fileInput.close()
 
 # Now generate the rest
-m = ones(numPart) * rho0 * boxSize**2 / numPart
+m = ones(numPart) * rho0 * boxSize ** 2 / numPart
 u = zeros(numPart)
 v = zeros((numPart, 3))
 
 for i in range(numPart):
-    
-    x = coords[i,0]
-    y = coords[i,1]
 
-    r2 = (x - boxSize / 2)**2 + (y - boxSize / 2)**2
+    x = coords[i, 0]
+    y = coords[i, 1]
+
+    r2 = (x - boxSize / 2) ** 2 + (y - boxSize / 2) ** 2
     r = sqrt(r2)
 
-    v_phi = 0.
+    v_phi = 0.0
     if r < 0.2:
-        v_phi = 5.*r
+        v_phi = 5.0 * r
     elif r < 0.4:
-        v_phi = 2. - 5.*r
+        v_phi = 2.0 - 5.0 * r
     else:
-        v_phi = 0.
-    v[i,0] = -v_phi * (y - boxSize / 2) / r
-    v[i,1] =  v_phi * (x - boxSize / 2) / r
-    v[i,2] = 0.
+        v_phi = 0.0
+    v[i, 0] = -v_phi * (y - boxSize / 2) / r
+    v[i, 1] = v_phi * (x - boxSize / 2) / r
+    v[i, 2] = 0.0
 
     P = P0
     if r < 0.2:
-        P = P + 5. + 12.5*r2
+        P = P + 5.0 + 12.5 * r2
     elif r < 0.4:
-        P = P + 9. + 12.5*r2 - 20.*r + 4.*log(r/0.2)
+        P = P + 9.0 + 12.5 * r2 - 20.0 * r + 4.0 * log(r / 0.2)
     else:
-        P = P + 3. + 4.*log(2.)
-    u[i] = P / ((gamma - 1.)*rho0)
-    
+        P = P + 3.0 + 4.0 * log(2.0)
+    u[i] = P / ((gamma - 1.0) * rho0)
 
 
-#File
-fileOutput = h5py.File(fileOutputName, 'w')
+# File
+fileOutput = h5py.File(fileOutputName, "w")
 
 # Header
 grp = fileOutput.create_group("/Header")
 grp.attrs["BoxSize"] = [boxSize, boxSize, 0.2]
-grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
 grp.attrs["Time"] = 0.0
@@ -89,29 +88,27 @@ grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["Dimension"] = 2
 
-#Units
+# Units
 grp = fileOutput.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 1.
-grp.attrs["Unit mass in cgs (U_M)"] = 1.
-grp.attrs["Unit time in cgs (U_t)"] = 1.
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+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
+# Particle group
 grp = fileOutput.create_group("/PartType0")
-ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
 ds[()] = coords
-ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ds = grp.create_dataset("Velocities", (numPart, 3), "f")
 ds[()] = v
-ds = grp.create_dataset('Masses', (numPart, 1), 'f')
-ds[()] = m.reshape((numPart,1))
-ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-ds[()] = h.reshape((numPart,1))
-ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-ds[()] = u.reshape((numPart,1))
-ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
-ds[()] = ids.reshape((numPart,1))
+ds = grp.create_dataset("Masses", (numPart, 1), "f")
+ds[()] = m.reshape((numPart, 1))
+ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+ds[()] = h.reshape((numPart, 1))
+ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+ds[()] = u.reshape((numPart, 1))
+ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+ds[()] = ids.reshape((numPart, 1))
 
 fileOutput.close()
-
-

examples/HydroTests/GreshoVortex_2D/plotSolution.py

 ###############################################################################
 # This file is part of SWIFT.
-# Copyright (c) 2016  Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+# Copyright (c) 2016  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
@@ -91,7 +91,7 @@ u = sim["/PartType0/InternalEnergies"][:]
 S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
 
-# Bin te data
+# Bin the data
 r_bin_edge = np.arange(0.0, 1.0, 0.02)
 r_bin = 0.5 * (r_bin_edge[1:] + r_bin_edge[:-1])
 rho_bin, _, _ = stats.binned_statistic(r, rho, statistic="mean", bins=r_bin_edge)

examples/HydroTests/GreshoVortex_2D/run.sh

@@ -9,11 +9,11 @@ fi
 if [ ! -e greshoVortex.hdf5 ]
 then
     echo "Generating initial conditions for the Gresho-Chan vortex example..."
-    python makeIC.py
+    python3 makeIC.py
 fi
 
 # Run SWIFT
-../../swift --hydro --threads=1 gresho.yml 2>&1 | tee output.log
+../../../swift --hydro --threads=1 gresho.yml 2>&1 | tee output.log
 
 # Plot the solution
-python plotSolution.py 11
+python3 plotSolution.py 11

examples/HydroTests/GreshoVortex_3D/getGlass.sh

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

examples/HydroTests/GreshoVortex_3D/makeIC.py

 ################################################################################
 # This file is part of SWIFT.
-# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+# Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
 #               2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
 #
 # This program is free software: you can redistribute it and/or modify
@@ -24,16 +24,16 @@ from numpy import *
 # Generates a swift IC file for the Gresho-Chan vortex in a periodic box
 
 # Parameters
-gamma = 5./3.     # Gas adiabatic index
-rho0 = 1           # Gas density
-P0 = 0.           # Constant additional pressure (should have no impact on the dynamics)
-fileOutputName = "greshoVortex.hdf5" 
+gamma = 5.0 / 3.0  # Gas adiabatic index
+rho0 = 1  # Gas density
+P0 = 0.0  # Constant additional pressure (should have no impact on the dynamics)
+fileOutputName = "greshoVortex.hdf5"
 fileGlass = "glassCube_64.hdf5"
-#---------------------------------------------------
+# ---------------------------------------------------
 
 # Get position and smoothing lengths from the glass
-fileInput = h5py.File(fileGlass, 'r')
-coords = fileInput["/PartType0/Coordinates"][:,:]
+fileInput = h5py.File(fileGlass, "r")
+coords = fileInput["/PartType0/Coordinates"][:, :]
 h = fileInput["/PartType0/SmoothingLength"][:]
 ids = fileInput["/PartType0/ParticleIDs"][:]
 boxSize = fileInput["/Header"].attrs["BoxSize"][0]
@@ -41,47 +41,46 @@ numPart = size(h)
 fileInput.close()
 
 # Now generate the rest
-m = ones(numPart) * rho0 * boxSize**3 / numPart
+m = ones(numPart) * rho0 * boxSize ** 3 / numPart
 u = zeros(numPart)
 v = zeros((numPart, 3))
 
 for i in range(numPart):
-    
-    x = coords[i,0]
-    y = coords[i,1]
 
-    r2 = (x - boxSize / 2)**2 + (y - boxSize / 2)**2
+    x = coords[i, 0]
+    y = coords[i, 1]
+
+    r2 = (x - boxSize / 2) ** 2 + (y - boxSize / 2) ** 2
     r = sqrt(r2)
 
-    v_phi = 0.
+    v_phi = 0.0
     if r < 0.2:
-        v_phi = 5.*r
+        v_phi = 5.0 * r
     elif r < 0.4:
-        v_phi = 2. - 5.*r
+        v_phi = 2.0 - 5.0 * r
     else:
-        v_phi = 0.
-    v[i,0] = -v_phi * (y - boxSize / 2) / r
-    v[i,1] =  v_phi * (x - boxSize / 2) / r
-    v[i,2] = 0.
+        v_phi = 0.0
+    v[i, 0] = -v_phi * (y - boxSize / 2) / r
+    v[i, 1] = v_phi * (x - boxSize / 2) / r
+    v[i, 2] = 0.0
 
     P = P0
     if r < 0.2:
-        P = P + 5. + 12.5*r2
+        P = P + 5.0 + 12.5 * r2
     elif r < 0.4:
-        P = P + 9. + 12.5*r2 - 20.*r + 4.*log(r/0.2)
+        P = P + 9.0 + 12.5 * r2 - 20.0 * r + 4.0 * log(r / 0.2)
     else:
-        P = P + 3. + 4.*log(2.)
-    u[i] = P / ((gamma - 1.)*rho0)
-    
+        P = P + 3.0 + 4.0 * log(2.0)
+    u[i] = P / ((gamma - 1.0) * rho0)
 
 
-#File
-fileOutput = h5py.File(fileOutputName, 'w')
+# File
+fileOutput = h5py.File(fileOutputName, "w")
 
 # Header
 grp = fileOutput.create_group("/Header")
 grp.attrs["BoxSize"] = [boxSize, boxSize, boxSize]
-grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
 grp.attrs["Time"] = 0.0
@@ -90,27 +89,27 @@ grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
 grp.attrs["Dimension"] = 3
 
-#Units
+# Units
 grp = fileOutput.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 1.
-grp.attrs["Unit mass in cgs (U_M)"] = 1.
-grp.attrs["Unit time in cgs (U_t)"] = 1.
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+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
+# Particle group
 grp = fileOutput.create_group("/PartType0")
-ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
 ds[()] = coords
-ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ds = grp.create_dataset("Velocities", (numPart, 3), "f")
 ds[()] = v
-ds = grp.create_dataset('Masses', (numPart, 1), 'f')
-ds[()] = m.reshape((numPart,1))
-ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-ds[()] = h.reshape((numPart,1))
-ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-ds[()] = u.reshape((numPart,1))
-ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
-ds[()] = ids.reshape((numPart,1))
+ds = grp.create_dataset("Masses", (numPart, 1), "f")
+ds[()] = m.reshape((numPart, 1))
+ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+ds[()] = h.reshape((numPart, 1))
+ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+ds[()] = u.reshape((numPart, 1))
+ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+ds[()] = ids.reshape((numPart, 1))
 
 fileOutput.close()

examples/HydroTests/GreshoVortex_3D/plotSolution.py

 ################################################################################
 # This file is part of SWIFT.
-# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+# Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
 #               2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
 #
 # This program is free software: you can redistribute it and/or modify
@@ -106,7 +106,7 @@ try:
 except:
     plot_viscosity = False
 
-# Bin te data
+# Bin the data
 r_bin_edge = np.arange(0.0, 1.0, 0.02)
 r_bin = 0.5 * (r_bin_edge[1:] + r_bin_edge[:-1])
 rho_bin, _, _ = stats.binned_statistic(r, rho, statistic="mean", bins=r_bin_edge)

examples/HydroTests/GreshoVortex_3D/run.sh

@@ -9,11 +9,11 @@ fi
 if [ ! -e greshoVortex.hdf5 ]
 then
     echo "Generating initial conditions for the Gresho-Chan vortex example..."
-    python makeIC.py
+    python3 makeIC.py
 fi
 
 # Run SWIFT
-../../swift --hydro --threads=4 gresho.yml 2>&1 | tee output.log
+../../../swift --hydro --threads=4 gresho.yml 2>&1 | tee output.log
 
 # Plot the solution
-python plotSolution.py 11
+python3 plotSolution.py 11

examples/HydroTests/InteractingBlastWaves_1D/getReference.sh

 #!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/interactingBlastWaves1D_exact.txt
+wget https://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/interactingBlastWaves1D_exact.txt