examples/HydroTests/InteractingBlastWaves_1D/makeIC.py

@@ -28,32 +28,32 @@ import h5py
 
 fileName = "interactingBlastWaves.hdf5"
 numPart = 800
-boxSize = 2.
+boxSize = 2.0
 
 coords = np.zeros((numPart, 3))
 v = np.zeros((numPart, 3))
 m = np.zeros(numPart) + boxSize / numPart
-h = np.zeros(numPart) + 2. * boxSize / numPart
+h = np.zeros(numPart) + 2.0 * boxSize / numPart
 u = np.zeros(numPart)
-ids = np.arange(numPart, dtype = 'L')
+ids = np.arange(numPart, dtype="L")
 rho = np.ones(numPart)
 
 for i in range(numPart):
-  coords[i,0] = (i + 0.5) * boxSize / numPart
-  if coords[i,0] < 0.1 or coords[i,0] > 1.9:
-    u[i] = 2500.
-  elif coords[i,0] > 0.9 and coords[i,0] < 1.1:
-    u[i] = 250.
-  else:
-    u[i] = 0.025
+    coords[i, 0] = (i + 0.5) * boxSize / numPart
+    if coords[i, 0] < 0.1 or coords[i, 0] > 1.9:
+        u[i] = 2500.0
+    elif coords[i, 0] > 0.9 and coords[i, 0] < 1.1:
+        u[i] = 250.0
+    else:
+        u[i] = 0.025
 
-#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"] =  [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
@@ -62,22 +62,22 @@ grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 grp.attrs["Flag_Entropy_ICs"] = 0
 grp.attrs["Dimension"] = 1
 
-#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=coords, 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('Density', data=rho, dtype='f')
+grp.create_dataset("Coordinates", data=coords, 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("Density", data=rho, dtype="f")
 
 file.close()

examples/HydroTests/InteractingBlastWaves_1D/plotSolution.py

@@ -17,116 +17,139 @@
 #
 ##############################################################################
 
-import numpy as np
 import matplotlib
+
 matplotlib.use("Agg")
-import pylab as pl
-import h5py
+import matplotlib.pyplot as plt
+import numpy as np
 import sys
+import h5py
+
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
-# Parameters
-gamma = 1.4 # Polytropic index
-
-# 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
-}
-pl.rcParams.update(params)
-pl.rc('font',**{'family':'sans-serif','sans-serif':['Times']})
-
-# Read the snapshot index from the command line argument
 snap = int(sys.argv[1])
 
 # Open the file and read the relevant data
 file = h5py.File("interactingBlastWaves_{0:04d}.hdf5".format(snap), "r")
-x = file["/PartType0/Coordinates"][:,0]
-rho = file["/PartType0/Densities"]
-v = file["/PartType0/Velocities"][:,0]
-u = file["/PartType0/InternalEnergies"]
-S = file["/PartType0/Entropies"]
-P = file["/PartType0/Pressures"]
+x = file["/PartType0/Coordinates"][:, 0]
+rho = file["/PartType0/Densities"][:]
+v = file["/PartType0/Velocities"][:, 0]
+u = file["/PartType0/InternalEnergies"][:]
+S = file["/PartType0/Entropies"][:]
+P = file["/PartType0/Pressures"][:]
 time = file["/Header"].attrs["Time"][0]
 
 scheme = file["/HydroScheme"].attrs["Scheme"]
 kernel = file["/HydroScheme"].attrs["Kernel function"]
 neighbours = file["/HydroScheme"].attrs["Kernel target N_ngb"][0]
 eta = file["/HydroScheme"].attrs["Kernel eta"][0]
+gamma = file["/HydroScheme"].attrs["Adiabatic index"]
 git = file["Code"].attrs["Git Revision"]
 
+if gamma != 1.4:
+    print(
+        "Error: SWIFT was run with the wrong adiabatic index. Should have been 1.4",
+        gamma,
+    )
+    exit(1)
+
 ref = np.loadtxt("interactingBlastWaves1D_exact.txt")
 
 # Plot the interesting quantities
-fig, ax = pl.subplots(2, 3)
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=4,
+    markeredgecolor="none",
+    alpha=0.2,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
 
 # Velocity profile
-ax[0][0].plot(x, v, "r.", markersize = 4.)
-ax[0][0].plot(ref[:,0], ref[:,2], "k--", alpha = 0.8, linewidth = 1.2)
-ax[0][0].set_xlabel("${\\rm{Position}}~x$", labelpad = 0)
-ax[0][0].set_ylabel("${\\rm{Velocity}}~v_x$", labelpad = 0)
-ax[0][0].set_xlim(0., 1.)
-ax[0][0].set_ylim(-1., 15.)
+plt.subplot(231)
+plt.plot(x, v, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 2], "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
+plt.xlim(0.0, 1.0)
+plt.ylim(-1.0, 15.0)
 
 # Density profile
-ax[0][1].plot(x, rho, "r.", markersize = 4.)
-ax[0][1].plot(ref[:,0], ref[:,1], "k--", alpha = 0.8, linewidth = 1.2)
-ax[0][1].set_xlabel("${\\rm{Position}}~x$", labelpad = 0)
-ax[0][1].set_ylabel("${\\rm{Density}}~\\rho$", labelpad = 0)
-ax[0][1].set_xlim(0., 1.)
+plt.subplot(232)
+plt.plot(x, rho, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 1], "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.xlim(0.0, 1.0)
 
 # Pressure profile
-ax[0][2].plot(x, P, "r.", markersize = 4.)
-ax[0][2].plot(ref[:,0], ref[:,3], "k--", alpha = 0.8, linewidth = 1.2)
-ax[0][2].set_xlabel("${\\rm{Position}}~x$", labelpad = 0)
-ax[0][2].set_ylabel("${\\rm{Pressure}}~P$", labelpad = 0)
-ax[0][2].set_xlim(0., 1.)
+plt.subplot(233)
+plt.plot(x, P, **scatter_props)
+plt.plot(ref[:, 0], ref[:, 3], "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.xlim(0.0, 1.0)
 
 # Internal energy profile
-ax[1][0].plot(x, u, "r.", markersize = 4.)
-ax[1][0].plot(ref[:,0], ref[:,3] / ref[:,1] / (gamma - 1.), "k--", alpha = 0.8,
-              linewidth = 1.2)
-ax[1][0].set_xlabel("${\\rm{Position}}~x$", labelpad = 0)
-ax[1][0].set_ylabel("${\\rm{Internal~Energy}}~u$", labelpad = 0)
-ax[1][0].set_xlim(0., 1.)
+plt.subplot(234)
+plt.plot(x, u, **scatter_props)
+plt.plot(
+    ref[:, 0],
+    ref[:, 3] / ref[:, 1] / (gamma - 1.0),
+    "--",
+    color=line_color,
+    alpha=0.8,
+    lw=1.2,
+)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.xlim(0.0, 1.0)
 
 # Entropy profile
-ax[1][1].plot(x, S, "r.", markersize = 4.)
-ax[1][1].plot(ref[:,0], ref[:,3] / ref[:,1]**gamma, "k--", alpha = 0.8,
-              linewidth = 1.2)
-ax[1][1].set_xlabel("${\\rm{Position}}~x$", labelpad = 0)
-ax[1][1].set_ylabel("${\\rm{Entropy}}~S$", labelpad = 0)
-ax[1][1].set_xlim(0., 1.)
+plt.subplot(235)
+plt.plot(x, S, **scatter_props)
+plt.plot(
+    ref[:, 0], ref[:, 3] / ref[:, 1] ** gamma, "--", color=line_color, alpha=0.8, lw=1.2
+)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+plt.xlim(0.0, 1.0)
 
 # Run information
-ax[1][2].set_frame_on(False)
-ax[1][2].text(-0.49, 0.9,
-  "Interacting blast waves test\nwith $\\gamma={0:.3f}$ in 1D at $t = {1:.2f}$".format(
-    gamma, time), fontsize = 10)
-ax[1][2].plot([-0.49, 0.1], [0.62, 0.62], "k-", lw = 1)
-ax[1][2].text(-0.49, 0.5, "$\\textsc{{Swift}}$ {0}".format(git), fontsize = 10)
-ax[1][2].text(-0.49, 0.4, scheme, fontsize = 10)
-ax[1][2].text(-0.49, 0.3, kernel, fontsize = 10)
-ax[1][2].text(-0.49, 0.2,
-  "${0:.2f}$ neighbours ($\\eta={1:.3f}$)".format(neighbours, eta),
-  fontsize = 10)
-ax[1][2].set_xlim(-0.5, 0.5)
-ax[1][2].set_ylim(0., 1.)
-ax[1][2].set_xticks([])
-ax[1][2].set_yticks([])
-
-pl.tight_layout()
-pl.savefig("InteractingBlastWaves.png", dpi = 200)
+plt.subplot(236, frameon=False)
+
+text_fontsize = 5
+
+plt.text(
+    -0.45,
+    0.9,
+    "Interacting blast waves test\nwith $\\gamma={0:.3f}$ in 1D at $t = {1:.2f}$".format(
+        gamma[0], 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.0, 1.0)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("InteractingBlastWaves.png", dpi=200)

examples/HydroTests/InteractingBlastWaves_1D/run.sh

@@ -4,11 +4,11 @@
 if [ ! -e interactingBlastWaves.hdf5 ]
 then
     echo "Generating initial conditions for the Sedov blast example..."
-    python makeIC.py
+    python3 makeIC.py
 fi
 
 # Run SWIFT
-../../swift --hydro --threads=1 interactingBlastWaves.yml 2>&1 | tee output.log
+../../../swift --hydro --threads=1  --limiter interactingBlastWaves.yml 2>&1 | tee output.log
 
 # Get the high resolution reference solution if not present.
 if [ ! -e interactingBlastWaves1D_exact.txt ]
@@ -18,4 +18,4 @@ then
 fi
 
 # Plot the solution
-python plotSolution.py 4
+python3 plotSolution.py 4

examples/HydroTests/KelvinHelmholtzGrowthRate_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/KelvinHelmholtzGrowthRate_2D/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,50 +24,56 @@ import numpy as np
 # Generates a swift IC file for the Kelvin-Helmholtz vortex in a periodic box
 
 # Parameters
-gamma = 5./3.     # Gas adiabatic index
-P0    = 2.5       # Pressure
-rho0  = 1.        # Density
-d     = 0.0317    # Thickness of the transition layer
-B     = 0.0005    # Amplitude of the seed velocity
+gamma = 5.0 / 3.0  # Gas adiabatic index
+P0 = 2.5  # Pressure
+rho0 = 1.0  # Density
+d = 0.0317  # Thickness of the transition layer
+B = 0.0005  # Amplitude of the seed velocity
 
 fileOutputName = "kelvinHelmholtzGrowthRate.hdf5"
 
-#---------------------------------------------------
+# ---------------------------------------------------
 
-glass = h5py.File("glassPlane_128.hdf5", 'r')
+glass = h5py.File("glassPlane_128.hdf5", "r")
 pos = glass["/PartType0/Coordinates"][:, :]
 h = glass["/PartType0/SmoothingLength"][:]
 
 N = len(h)
-vol = 1.
+vol = 1.0
 
 # Generate extra arrays
 v = np.zeros((N, 3))
 ids = np.linspace(1, N, N)
 m = np.ones(N) * rho0 * vol / N
-u = np.ones(N) * P0 / (rho0 * (gamma - 1.))
+u = np.ones(N) * P0 / (rho0 * (gamma - 1.0))
 
 v[pos[:, 1] <= 0.25 - d, 0] = -0.5
-v[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 0] = \
-  -0.5 + \
-  0.5 * (pos[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 1] + d - 0.25) / d
+v[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 0] = (
+    -0.5
+    + 0.5 * (pos[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 1] + d - 0.25) / d
+)
 v[(pos[:, 1] <= 0.75 - d) & (pos[:, 1] >= 0.25 + d), 0] = 0.5
-v[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 0] = \
-  0.5 - \
-  0.5 * (pos[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 1] + d - 0.75) / d
+v[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 0] = (
+    0.5 - 0.5 * (pos[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 1] + d - 0.75) / d
+)
 v[pos[:, 1] >= 0.75 + d, 0] = -0.5
 
-v[:, 1] = B * np.sin(4. * np.pi * pos[:, 0]) * \
-          (np.exp(-(pos[:, 1] - 0.25)**2 / 32. / d**2) + \
-           np.exp(-(pos[:, 1] - 0.75)**2 / 32. / d**2))
-            
-#File
-fileOutput = h5py.File(fileOutputName, 'w')
+v[:, 1] = (
+    B
+    * np.sin(4.0 * np.pi * pos[:, 0])
+    * (
+        np.exp(-(pos[:, 1] - 0.25) ** 2 / 32.0 / d ** 2)
+        + np.exp(-(pos[:, 1] - 0.75) ** 2 / 32.0 / d ** 2)
+    )
+)
+
+# File
+fileOutput = h5py.File(fileOutputName, "w")
 
 # Header
 grp = fileOutput.create_group("/Header")
-grp.attrs["BoxSize"] = [1., 1., 1.]
-grp.attrs["NumPart_Total"] =  [N, 0, 0, 0, 0, 0]
+grp.attrs["BoxSize"] = [1.0, 1.0, 1.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
@@ -76,21 +82,21 @@ 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")
-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")
 
 fileOutput.close()

examples/HydroTests/KelvinHelmholtzGrowthRate_2D/makeIC_regular.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,56 +24,62 @@ import numpy as np
 # Generates a swift IC file for the Kelvin-Helmholtz vortex in a periodic box
 
 # Parameters
-gamma = 5./3.     # Gas adiabatic index
-P0    = 2.5       # Pressure
-rho0  = 1.        # Density
-d     = 0.0317    # Thickness of the transition layer
-B     = 0.0005    # Amplitude of the seed velocity
-N1D   = 128       # Number of particles in one dimension
+gamma = 5.0 / 3.0  # Gas adiabatic index
+P0 = 2.5  # Pressure
+rho0 = 1.0  # Density
+d = 0.0317  # Thickness of the transition layer
+B = 0.0005  # Amplitude of the seed velocity
+N1D = 128  # Number of particles in one dimension
 
 fileOutputName = "kelvinHelmholtzGrowthRate.hdf5"
 
-#---------------------------------------------------
+# ---------------------------------------------------
 
 N = N1D ** 2
-x = np.linspace(0., 1., N1D + 1)
+x = np.linspace(0.0, 1.0, N1D + 1)
 x = 0.5 * (x[1:] + x[:-1])
 y = x
 xx, yy = np.meshgrid(x, y)
 pos = np.zeros((N, 3))
 pos[:, 0] = xx.reshape((N))
 pos[:, 1] = yy.reshape((N))
-h = np.ones(N) * 2. / N1D
+h = np.ones(N) * 2.0 / N1D
 
-vol = 1.
+vol = 1.0
 
 # Generate extra arrays
 v = np.zeros((N, 3))
 ids = np.linspace(1, N, N)
 m = np.ones(N) * rho0 * vol / N
-u = np.ones(N) * P0 / (rho0 * (gamma - 1.))
+u = np.ones(N) * P0 / (rho0 * (gamma - 1.0))
 
 v[pos[:, 1] <= 0.25 - d, 0] = -0.5
-v[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 0] = \
-  -0.5 + \
-  0.5 * (pos[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 1] + d - 0.25) / d
+v[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 0] = (
+    -0.5
+    + 0.5 * (pos[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 1] + d - 0.25) / d
+)
 v[(pos[:, 1] <= 0.75 - d) & (pos[:, 1] >= 0.25 + d), 0] = 0.5
-v[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 0] = \
-  0.5 - \
-  0.5 * (pos[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 1] + d - 0.75) / d
+v[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 0] = (
+    0.5 - 0.5 * (pos[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 1] + d - 0.75) / d
+)
 v[pos[:, 1] >= 0.75 + d, 0] = -0.5
 
-v[:, 1] = B * np.sin(4. * np.pi * pos[:, 0]) * \
-          (np.exp(-(pos[:, 1] - 0.25)**2 / 32. / d**2) + \
-           np.exp(-(pos[:, 1] - 0.75)**2 / 32. / d**2))
-            
-#File
-fileOutput = h5py.File(fileOutputName, 'w')
+v[:, 1] = (
+    B
+    * np.sin(4.0 * np.pi * pos[:, 0])
+    * (
+        np.exp(-(pos[:, 1] - 0.25) ** 2 / 32.0 / d ** 2)
+        + np.exp(-(pos[:, 1] - 0.75) ** 2 / 32.0 / d ** 2)
+    )
+)
+
+# File
+fileOutput = h5py.File(fileOutputName, "w")
 
 # Header
 grp = fileOutput.create_group("/Header")
-grp.attrs["BoxSize"] = [1., 1., 1.]
-grp.attrs["NumPart_Total"] =  [N, 0, 0, 0, 0, 0]
+grp.attrs["BoxSize"] = [1.0, 1.0, 1.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
@@ -82,21 +88,21 @@ 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")
-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")
 
 fileOutput.close()

examples/HydroTests/KelvinHelmholtzGrowthRate_2D/plotSolution.py

@@ -20,28 +20,29 @@
 import numpy as np
 import h5py
 import matplotlib
+
 matplotlib.use("Agg")
 import pylab as pl
 import sys
 
 if len(sys.argv) < 2:
-  print "No final snapshot number provided!"
-  exit()
+    print("No final snapshot number provided!")
+    exit()
 lastsnap = int(sys.argv[1])
 
 # Read the simulation data
 t = np.zeros(lastsnap + 1)
 ey = np.zeros(lastsnap + 1)
 for i in range(lastsnap + 1):
-  file = h5py.File("kelvinHelmholtzGrowthRate_{0:04d}.hdf5".format(i), 'r')
-  t_snap = float(file["/Header"].attrs["Time"])
-  vy = file["/PartType0/Velocities"][:,1]
-  m = file["/PartType0/Masses"][:]
+    file = h5py.File("kelvinHelmholtzGrowthRate_{0:04d}.hdf5".format(i), "r")
+    t_snap = float(file["/Header"].attrs["Time"])
+    vy = file["/PartType0/Velocities"][:, 1]
+    m = file["/PartType0/Masses"][:]
 
-  ey_snap = 0.5 * m * vy**2
+    ey_snap = 0.5 * m * vy ** 2
 
-  t[i] = t_snap
-  ey[i] = ey_snap.sum()
+    t[i] = t_snap
+    ey[i] = ey_snap.sum()
 
 pl.semilogy(t, ey, "k.")
 pl.savefig("kelvinHelmholtzGrowthRate.png")

examples/HydroTests/KelvinHelmholtzGrowthRate_2D/run.sh

 #!/bin/bash
 
- # Generate the initial conditions if they are not present.
+# Generate the initial conditions if they are not present.
+if [ ! -e glassPlane_128.hdf5 ]
+then
+    ./getGlass.sh
+fi
+
 if [ ! -e kelvinHelmholtzGrowthRate.hdf5 ]
 then
     echo "Generating initial conditions for the Kelvin-Helmholtz growth rate " \
          "example..."
-    python makeIC.py
+    python3 makeIC.py
 fi
 
 # Run SWIFT
-../../swift --hydro --threads=1 kelvinHelmholtzGrowthRate.yml 2>&1 | tee output.log
+../../../swift --hydro --threads=1 kelvinHelmholtzGrowthRate.yml 2>&1 | tee output.log
 
 # Plot the solution
-python plotSolution.py 100
+python3 plotSolution.py 100

examples/HydroTests/KelvinHelmholtzGrowthRate_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/KelvinHelmholtzGrowthRate_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,50 +24,56 @@ import numpy as np
 # Generates a swift IC file for the Kelvin-Helmholtz vortex in a periodic box
 
 # Parameters
-gamma = 5./3.     # Gas adiabatic index
-P0    = 2.5       # Pressure
-rho0  = 1.        # Density
-d     = 0.0317    # Thickness of the transition layer
-B     = 0.0005    # Amplitude of the seed velocity
+gamma = 5.0 / 3.0  # Gas adiabatic index
+P0 = 2.5  # Pressure
+rho0 = 1.0  # Density
+d = 0.0317  # Thickness of the transition layer
+B = 0.0005  # Amplitude of the seed velocity
 
 fileOutputName = "kelvinHelmholtzGrowthRate.hdf5"
 
-#---------------------------------------------------
+# ---------------------------------------------------
 
-glass = h5py.File("glassCube_64.hdf5", 'r')
+glass = h5py.File("glassCube_64.hdf5", "r")
 pos = glass["/PartType0/Coordinates"][:, :]
 h = glass["/PartType0/SmoothingLength"][:]
 
 N = len(h)
-vol = 1.
+vol = 1.0
 
 # Generate extra arrays
 v = np.zeros((N, 3))
 ids = np.linspace(1, N, N)
 m = np.ones(N) * rho0 * vol / N
-u = np.ones(N) * P0 / (rho0 * (gamma - 1.))
+u = np.ones(N) * P0 / (rho0 * (gamma - 1.0))
 
 v[pos[:, 1] <= 0.25 - d, 0] = -0.5
-v[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 0] = \
-  -0.5 + \
-  0.5 * (pos[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 1] + d - 0.25) / d
+v[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 0] = (
+    -0.5
+    + 0.5 * (pos[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 1] + d - 0.25) / d
+)
 v[(pos[:, 1] <= 0.75 - d) & (pos[:, 1] >= 0.25 + d), 0] = 0.5
-v[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 0] = \
-  0.5 - \
-  0.5 * (pos[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 1] + d - 0.75) / d
+v[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 0] = (
+    0.5 - 0.5 * (pos[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 1] + d - 0.75) / d
+)
 v[pos[:, 1] >= 0.75 + d, 0] = -0.5
 
-v[:, 1] = B * np.sin(4. * np.pi * pos[:, 0]) * \
-          (np.exp(-(pos[:, 1] - 0.25)**2 / 32. / d**2) + \
-           np.exp(-(pos[:, 1] - 0.75)**2 / 32. / d**2))
-            
-#File
-fileOutput = h5py.File(fileOutputName, 'w')
+v[:, 1] = (
+    B
+    * np.sin(4.0 * np.pi * pos[:, 0])
+    * (
+        np.exp(-(pos[:, 1] - 0.25) ** 2 / 32.0 / d ** 2)
+        + np.exp(-(pos[:, 1] - 0.75) ** 2 / 32.0 / d ** 2)
+    )
+)
+
+# File
+fileOutput = h5py.File(fileOutputName, "w")
 
 # Header
 grp = fileOutput.create_group("/Header")
-grp.attrs["BoxSize"] = [1., 1., 1.]
-grp.attrs["NumPart_Total"] =  [N, 0, 0, 0, 0, 0]
+grp.attrs["BoxSize"] = [1.0, 1.0, 1.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
@@ -76,21 +82,21 @@ 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")
-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")
 
 fileOutput.close()

examples/HydroTests/KelvinHelmholtzGrowthRate_3D/makeIC_regular.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,19 +24,19 @@ import numpy as np
 # Generates a swift IC file for the Kelvin-Helmholtz vortex in a periodic box
 
 # Parameters
-gamma = 5./3.     # Gas adiabatic index
-P0    = 2.5       # Pressure
-rho0  = 1.        # Density
-d     = 0.0317    # Thickness of the transition layer
-B     = 0.0005    # Amplitude of the seed velocity
-N1D   = 64        # Number of particles in one dimension
+gamma = 5.0 / 3.0  # Gas adiabatic index
+P0 = 2.5  # Pressure
+rho0 = 1.0  # Density
+d = 0.0317  # Thickness of the transition layer
+B = 0.0005  # Amplitude of the seed velocity
+N1D = 64  # Number of particles in one dimension
 
 fileOutputName = "kelvinHelmholtzGrowthRate.hdf5"
 
-#---------------------------------------------------
+# ---------------------------------------------------
 
 N = N1D ** 3
-x = np.linspace(0., 1., N1D + 1)
+x = np.linspace(0.0, 1.0, N1D + 1)
 x = 0.5 * (x[1:] + x[:-1])
 y = x
 z = x
@@ -45,37 +45,43 @@ pos = np.zeros((N, 3))
 pos[:, 0] = xx.reshape((N))
 pos[:, 1] = yy.reshape((N))
 pos[:, 2] = zz.reshape((N))
-h = np.ones(N) * 2. / N1D
+h = np.ones(N) * 2.0 / N1D
 
-vol = 1.
+vol = 1.0
 
 # Generate extra arrays
 v = np.zeros((N, 3))
 ids = np.linspace(1, N, N)
 m = np.ones(N) * rho0 * vol / N
-u = np.ones(N) * P0 / (rho0 * (gamma - 1.))
+u = np.ones(N) * P0 / (rho0 * (gamma - 1.0))
 
 v[pos[:, 1] <= 0.25 - d, 0] = -0.5
-v[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 0] = \
-  -0.5 + \
-  0.5 * (pos[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 1] + d - 0.25) / d
+v[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 0] = (
+    -0.5
+    + 0.5 * (pos[(pos[:, 1] < 0.25 + d) & (pos[:, 1] > 0.25 - d), 1] + d - 0.25) / d
+)
 v[(pos[:, 1] <= 0.75 - d) & (pos[:, 1] >= 0.25 + d), 0] = 0.5
-v[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 0] = \
-  0.5 - \
-  0.5 * (pos[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 1] + d - 0.75) / d
+v[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 0] = (
+    0.5 - 0.5 * (pos[(pos[:, 1] < 0.75 + d) & (pos[:, 1] > 0.75 - d), 1] + d - 0.75) / d
+)
 v[pos[:, 1] >= 0.75 + d, 0] = -0.5
 
-v[:, 1] = B * np.sin(4. * np.pi * pos[:, 0]) * \
-          (np.exp(-(pos[:, 1] - 0.25)**2 / 32. / d**2) + \
-           np.exp(-(pos[:, 1] - 0.75)**2 / 32. / d**2))
-            
-#File
-fileOutput = h5py.File(fileOutputName, 'w')
+v[:, 1] = (
+    B
+    * np.sin(4.0 * np.pi * pos[:, 0])
+    * (
+        np.exp(-(pos[:, 1] - 0.25) ** 2 / 32.0 / d ** 2)
+        + np.exp(-(pos[:, 1] - 0.75) ** 2 / 32.0 / d ** 2)
+    )
+)
+
+# File
+fileOutput = h5py.File(fileOutputName, "w")
 
 # Header
 grp = fileOutput.create_group("/Header")
-grp.attrs["BoxSize"] = [1., 1., 1.]
-grp.attrs["NumPart_Total"] =  [N, 0, 0, 0, 0, 0]
+grp.attrs["BoxSize"] = [1.0, 1.0, 1.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
@@ -84,21 +90,21 @@ 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")
-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")
 
 fileOutput.close()

examples/HydroTests/KelvinHelmholtzGrowthRate_3D/plotSolution.py

@@ -20,28 +20,29 @@
 import numpy as np
 import h5py
 import matplotlib
+
 matplotlib.use("Agg")
 import pylab as pl
 import sys
 
 if len(sys.argv) < 2:
-  print "No final snapshot number provided!"
-  exit()
+    print("No final snapshot number provided!")
+    exit()
 lastsnap = int(sys.argv[1])
 
 # Read the simulation data
 t = np.zeros(lastsnap + 1)
 ey = np.zeros(lastsnap + 1)
 for i in range(lastsnap + 1):
-  file = h5py.File("kelvinHelmholtzGrowthRate_{0:04d}.hdf5".format(i), 'r')
-  t_snap = float(file["/Header"].attrs["Time"])
-  vy = file["/PartType0/Velocities"][:,1]
-  m = file["/PartType0/Masses"][:]
+    file = h5py.File("kelvinHelmholtzGrowthRate_{0:04d}.hdf5".format(i), "r")
+    t_snap = float(file["/Header"].attrs["Time"])
+    vy = file["/PartType0/Velocities"][:, 1]
+    m = file["/PartType0/Masses"][:]
 
-  ey_snap = 0.5 * m * vy**2
+    ey_snap = 0.5 * m * vy ** 2
 
-  t[i] = t_snap
-  ey[i] = ey_snap.sum()
+    t[i] = t_snap
+    ey[i] = ey_snap.sum()
 
 pl.semilogy(t, ey, "k.")
 pl.savefig("kelvinHelmholtzGrowthRate.png")

examples/HydroTests/KelvinHelmholtzGrowthRate_3D/run.sh

 #!/bin/bash
 
- # Generate the initial conditions if they are not present.
+# Generate the initial conditions if they are not present.
+if [ ! -e glassCube_64.hdf5 ]
+then
+    ./getGlass.sh
+fi
+
 if [ ! -e kelvinHelmholtzGrowthRate.hdf5 ]
 then
     echo "Generating initial conditions for the Kelvin-Helmholtz growth rate " \
          "example..."
-    python makeIC.py
+    python3 makeIC.py
 fi
 
 # Run SWIFT
-../../swift --hydro --threads=1 kelvinHelmholtzGrowthRate.yml 2>&1 | tee output.log
+../../../swift --hydro --threads=1 kelvinHelmholtzGrowthRate.yml 2>&1 | tee output.log
 
 # Plot the solution
-python plotSolution.py 100
+python3 plotSolution.py 100

examples/HydroTests/KelvinHelmholtz_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 *
@@ -24,27 +24,27 @@ import sys
 # Generates a swift IC file for the Kelvin-Helmholtz vortex in a periodic box
 
 # Parameters
-L2    = 256       # Particles along one edge in the low-density region
-gamma = 5./3.     # Gas adiabatic index
-P1    = 2.5       # Central region pressure
-P2    = 2.5       # Outskirts pressure
-v1    = 0.5       # Central region velocity
-v2    = -0.5      # Outskirts vlocity
-rho1  = 2         # Central density
-rho2  = 1         # Outskirts density
+L2 = 256  # Particles along one edge in the low-density region
+gamma = 5.0 / 3.0  # Gas adiabatic index
+P1 = 2.5  # Central region pressure
+P2 = 2.5  # Outskirts pressure
+v1 = 0.5  # Central region velocity
+v2 = -0.5  # Outskirts vlocity
+rho1 = 2  # Central density
+rho2 = 1  # Outskirts density
 omega0 = 0.1
 sigma = 0.05 / sqrt(2)
-fileOutputName = "kelvinHelmholtz.hdf5" 
-#---------------------------------------------------
+fileOutputName = "kelvinHelmholtz.hdf5"
+# ---------------------------------------------------
 
 # Start by generating grids of particles at the two densities
 numPart2 = L2 * L2
 L1 = int(sqrt(numPart2 / rho2 * rho1))
 numPart1 = L1 * L1
 
-#print "N2 =", numPart2, "N1 =", numPart1
-#print "L2 =", L2, "L1 = ", L1
-#print "rho2 =", rho2, "rho1 =", (float(L1*L1)) / (float(L2*L2))
+print("N2 =", numPart2, "N1 =", numPart1)
+print("L2 =", L2, "L1 = ", L1)
+print("rho2 =", rho2, "rho1 =", (float(L1 * L1)) / (float(L2 * L2)))
 
 coords1 = zeros((numPart1, 3))
 coords2 = zeros((numPart2, 3))
@@ -62,39 +62,39 @@ for i in range(L1):
     for j in range(L1):
 
         index = i * L1 + j
-        
-        x = i / float(L1) + 1. / (2. * L1)
-        y = j / float(L1) + 1. / (2. * L1)
+
+        x = i / float(L1) + 1.0 / (2.0 * L1)
+        y = j / float(L1) + 1.0 / (2.0 * L1)
 
         coords1[index, 0] = x
         coords1[index, 1] = y
-        u1[index] = P1 / (rho1 * (gamma-1.))
+        u1[index] = P1 / (rho1 * (gamma - 1.0))
         vel1[index, 0] = v1
-        
+
 # Particles in the outskirts
 for i in range(L2):
     for j in range(L2):
 
         index = i * L2 + j
-        
-        x = i / float(L2) + 1. / (2. * L2)
-        y = j / float(L2) + 1. / (2. * L2)
+
+        x = i / float(L2) + 1.0 / (2.0 * L2)
+        y = j / float(L2) + 1.0 / (2.0 * L2)
 
         coords2[index, 0] = x
         coords2[index, 1] = y
-        u2[index] = P2 / (rho2 * (gamma-1.))
+        u2[index] = P2 / (rho2 * (gamma - 1.0))
         vel2[index, 0] = v2
 
 
 # Now concatenate arrays
-where1 = abs(coords1[:,1]-0.5) < 0.25
-where2 = abs(coords2[:,1]-0.5) > 0.25
+where1 = abs(coords1[:, 1] - 0.5) < 0.25
+where2 = abs(coords2[:, 1] - 0.5) > 0.25
 
 coords = append(coords1[where1, :], coords2[where2, :], axis=0)
 
-#print L2*(L2/2), L1*(L1/2)
-#print shape(coords), shape(coords1[where1,:]), shape(coords2[where2,:])
-#print shape(coords), shape(logical_not(coords1[where1,:])), shape(logical_not(coords2[where2,:]))
+# print L2*(L2/2), L1*(L1/2)
+# print shape(coords), shape(coords1[where1,:]), shape(coords2[where2,:])
+# print shape(coords), shape(logical_not(coords1[where1,:])), shape(logical_not(coords2[where2,:]))
 
 vel = append(vel1[where1, :], vel2[where2, :], axis=0)
 h = append(h1[where1], h2[where2], axis=0)
@@ -105,15 +105,22 @@ ids = linspace(1, numPart, numPart)
 m[:] = (0.5 * rho1 + 0.5 * rho2) / float(numPart)
 
 # Velocity perturbation
-vel[:,1] = omega0 * sin(4*pi*coords[:,0]) * (exp(-(coords[:,1]-0.25)**2 / (2 * sigma**2)) + exp(-(coords[:,1]-0.75)**2 / (2 * sigma**2)))
-            
-#File
-fileOutput = h5py.File(fileOutputName, 'w')
+vel[:, 1] = (
+    omega0
+    * sin(4 * pi * coords[:, 0])
+    * (
+        exp(-((coords[:, 1] - 0.25) ** 2) / (2 * sigma ** 2))
+        + exp(-((coords[:, 1] - 0.75) ** 2) / (2 * sigma ** 2))
+    )
+)
+
+# File
+fileOutput = h5py.File(fileOutputName, "w")
 
 # Header
 grp = fileOutput.create_group("/Header")
-grp.attrs["BoxSize"] = [1., 1., 0.1]
-grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["BoxSize"] = [1.0, 1.0, 0.1]
+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
@@ -122,29 +129,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[()] = vel
-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/KelvinHelmholtz_2D/makeMovie.py

@@ -46,7 +46,7 @@ def make_plot(filename, array, nx, ny, dx, dy):
 
     array.set_array(mesh)
 
-    return array,
+    return (array,)
 
 
 def frame(n, *args):
@@ -63,6 +63,7 @@ def frame(n, *args):
 
 if __name__ == "__main__":
     import matplotlib
+
     matplotlib.use("Agg")
 
     from tqdm import tqdm
@@ -74,7 +75,6 @@ if __name__ == "__main__":
     filename = "kelvinHelmholtz"
     dpi = 512
 
-
     # Look for the number of files in the directory.
     i = 0
     while True:
@@ -84,8 +84,7 @@ if __name__ == "__main__":
             break
 
         if i > 10000:
-            raise FileNotFoundError(
-                "Could not find the snapshots in the directory")
+            raise FileNotFoundError("Could not find the snapshots in the directory")
 
     frames = tqdm(np.arange(0, i))
 
@@ -100,7 +99,7 @@ if __name__ == "__main__":
     xv, yv = np.meshgrid(x, y)
 
     mesh = si.griddata((data_x, data_y), density, (xv, yv), method="nearest")
-    
+
     # Global variable for set_array
     plot = ax.imshow(mesh, extent=[0, 1, 0, 1], animated=True, interpolation="none")
 

examples/HydroTests/KelvinHelmholtz_2D/plotSolution.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/>.
+#
+##############################################################################
 
 # Computes the analytical solution of the Gresho-Chan vortex and plots the SPH answer
 
 # Parameters
-gas_gamma = 5./3.     # Gas adiabatic index
-P1    = 2.5       # Central region pressure
-P2    = 2.5       # Outskirts pressure
-v1    = 0.5       # Central region velocity
-v2    = -0.5      # Outskirts vlocity
-rho1  = 2         # Central density
-rho2  = 1         # Outskirts density
+gas_gamma = 5.0 / 3.0  # Gas adiabatic index
+P1 = 2.5  # Central region pressure
+P2 = 2.5  # Outskirts pressure
+v1 = 0.5  # Central region velocity
+v2 = -0.5  # Outskirts vlocity
+rho1 = 2  # Central density
+rho2 = 1  # Outskirts density
 
 # ---------------------------------------------------------------
 # Don't touch anything after this.
 # ---------------------------------------------------------------
 
 import matplotlib
+
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 import h5py
+import sys
 
-# 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']})
-
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
 # Read the simulation data
-sim = h5py.File("kelvinHelmholtz_%04d.hdf5"%snap, "r")
+sim = h5py.File("kelvinHelmholtz_%04d.hdf5" % snap, "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
 time = sim["/Header"].attrs["Time"][0]
 scheme = sim["/HydroScheme"].attrs["Scheme"]
@@ -72,88 +54,135 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
 eta = sim["/HydroScheme"].attrs["Kernel eta"]
 git = sim["Code"].attrs["Git Revision"]
 
-pos = sim["/PartType0/Coordinates"][:,:]
-x = pos[:,0] - boxSize / 2
-y = pos[:,1] - boxSize / 2
-vel = sim["/PartType0/Velocities"][:,:]
-v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
+pos = sim["/PartType0/Coordinates"][:, :]
+x = pos[:, 0] - boxSize / 2
+y = pos[:, 1] - boxSize / 2
+vel = sim["/PartType0/Velocities"][:, :]
+v_norm = np.sqrt(vel[:, 0] ** 2 + vel[:, 1] ** 2)
 rho = sim["/PartType0/Densities"][:]
 u = sim["/PartType0/InternalEnergies"][:]
 S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
 
 # Plot the interesting quantities
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
 
 
 # Azimuthal velocity profile -----------------------------
-subplot(231)
-scatter(pos[:,0], pos[:,1], c=vel[:,0], cmap="PuBu", edgecolors='face', s=4, vmin=-1., vmax=1.)
-text(0.97, 0.97, "${\\rm{Velocity~along}}~x$", ha="right", va="top", backgroundcolor="w")
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=0)
-xlim(0, 1)
-ylim(0, 1)
+plt.subplot(231)
+plt.scatter(
+    pos[:, 0],
+    pos[:, 1],
+    c=vel[:, 0],
+    cmap="PuBu",
+    edgecolors="face",
+    s=0.25,
+    vmin=-1.0,
+    vmax=1.0,
+)
+plt.text(
+    0.97, 0.97, "${\\rm{Velocity~along}}~x$", ha="right", va="top", backgroundcolor="w"
+)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Position}}~y$", labelpad=0)
+plt.xlim(0, 1)
+plt.ylim(0, 1)
 
 # Radial density profile --------------------------------
-subplot(232)
-scatter(pos[:,0], pos[:,1], c=rho, cmap="PuBu", edgecolors='face', s=4, vmin=0.8, vmax=2.2)
-text(0.97, 0.97, "${\\rm{Density}}$", ha="right", va="top", backgroundcolor="w")
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=0)
-xlim(0, 1)
-ylim(0, 1)
+plt.subplot(232)
+plt.scatter(
+    pos[:, 0],
+    pos[:, 1],
+    c=rho,
+    cmap="PuBu",
+    edgecolors="face",
+    s=0.25,
+    vmin=0.8,
+    vmax=2.2,
+)
+plt.text(0.97, 0.97, "${\\rm{Density}}$", ha="right", va="top", backgroundcolor="w")
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Position}}~y$", labelpad=0)
+plt.xlim(0, 1)
+plt.ylim(0, 1)
 
 # Radial pressure profile --------------------------------
-subplot(233)
-scatter(pos[:,0], pos[:,1], c=P, cmap="PuBu", edgecolors='face', s=4, vmin=1, vmax=4)
-text(0.97, 0.97, "${\\rm{Pressure}}$", ha="right", va="top", backgroundcolor="w")
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=0)
-xlim(0, 1)
-ylim(0, 1)
+plt.subplot(233)
+plt.scatter(
+    pos[:, 0], pos[:, 1], c=P, cmap="PuBu", edgecolors="face", s=0.25, vmin=1, vmax=4
+)
+plt.text(0.97, 0.97, "${\\rm{Pressure}}$", ha="right", va="top", backgroundcolor="w")
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Position}}~y$", labelpad=0)
+plt.xlim(0, 1)
+plt.ylim(0, 1)
 
 # Internal energy profile --------------------------------
-subplot(234)
-scatter(pos[:,0], pos[:,1], c=u, cmap="PuBu", edgecolors='face', s=4, vmin=1.5, vmax=5.)
-text(0.97, 0.97, "${\\rm{Internal~energy}}$", ha="right", va="top", backgroundcolor="w")
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=0)
-xlim(0, 1)
-ylim(0, 1)
+plt.subplot(234)
+plt.scatter(
+    pos[:, 0],
+    pos[:, 1],
+    c=u,
+    cmap="PuBu",
+    edgecolors="face",
+    s=0.25,
+    vmin=1.5,
+    vmax=5.0,
+)
+plt.text(
+    0.97, 0.97, "${\\rm{Internal~energy}}$", ha="right", va="top", backgroundcolor="w"
+)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Position}}~y$", labelpad=0)
+plt.xlim(0, 1)
+plt.ylim(0, 1)
 
 # Radial entropy profile --------------------------------
-subplot(235)
-scatter(pos[:,0], pos[:,1], c=S, cmap="PuBu", edgecolors='face', s=4, vmin=0.5, vmax=3.)
-text(0.97, 0.97, "${\\rm{Entropy}}$", ha="right", va="top", backgroundcolor="w")
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=0)
-xlim(0, 1)
-ylim(0, 1)
-
-# Image --------------------------------------------------
-#subplot(234)
-#scatter(pos[:,0], pos[:,1], c=v_norm, cmap="PuBu", edgecolors='face', s=4, vmin=0, vmax=1)
-#text(0.95, 0.95, "$|v|$", ha="right", va="top")
-#xlim(0,1)
-#ylim(0,1)
-#xlabel("$x$", labelpad=0)
-#ylabel("$y$", labelpad=0)
+plt.subplot(235)
+plt.scatter(
+    pos[:, 0],
+    pos[:, 1],
+    c=S,
+    cmap="PuBu",
+    edgecolors="face",
+    s=0.25,
+    vmin=0.5,
+    vmax=3.0,
+)
+plt.text(0.97, 0.97, "${\\rm{Entropy}}$", ha="right", va="top", backgroundcolor="w")
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Position}}~y$", labelpad=0)
+plt.xlim(0, 1)
+plt.ylim(0, 1)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
-
-text(-0.49, 0.9, "Kelvin-Helmholtz instability at $t=%.2f$"%(time), fontsize=10)
-text(-0.49, 0.8, "Centre:~~~ $(P, \\rho, v) = (%.3f, %.3f, %.3f)$"%(P1, rho1, v1), fontsize=10)
-text(-0.49, 0.7, "Outskirts: $(P, \\rho, v) = (%.3f, %.3f, %.3f)$"%(P2, rho2, v2), 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([])
-
-savefig("KelvinHelmholtz.png", dpi=200)
+plt.subplot(236, frameon=False)
+
+plt.text(-0.45, 0.9, "Kelvin-Helmholtz instability at $t=%.2f$" % (time), fontsize=10)
+plt.text(
+    -0.45,
+    0.8,
+    "Centre: $(P, \\rho, v) = (%.3f, %.3f, %.3f)$" % (P1, rho1, v1),
+    fontsize=10,
+)
+plt.text(
+    -0.45,
+    0.7,
+    "Outskirts: $(P, \\rho, v) = (%.3f, %.3f, %.3f)$" % (P2, rho2, v2),
+    fontsize=10,
+)
+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=10)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=10)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=10)
+plt.text(
+    -0.45, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta), fontsize=10
+)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("KelvinHelmholtz.png", dpi=200)

examples/HydroTests/KelvinHelmholtz_2D/run.sh

@@ -4,12 +4,13 @@
 if [ ! -e kelvinHelmholtz.hdf5 ]
 then
     echo "Generating initial conditions for the Kelvin-Helmholtz example..."
-    python makeIC.py
+    python3 makeIC.py
 fi
 
 # Run SWIFT
-../../swift --hydro --threads=4 kelvinHelmholtz.yml 2>&1 | tee output.log
+../../../swift --hydro --threads=4 kelvinHelmholtz.yml 2>&1 | tee output.log
 
 
 # Plot the solution
+python3 plotSolution.py 100
 python3 makeMovieSwiftsimIO.py

examples/HydroTests/KeplerianRing/README.md

@@ -18,9 +18,6 @@ The test uses:
 + $c_s = 0.01 \ll v_\phi = 10$, enforced by giving them a mass of 1 unit and an
   internal energy of 0.015.
 
-Please note that the initial condition generator **requires python3 rather than
-python2**, as well as the plotting script.
-
 
 Code Setup
 ----------
@@ -69,7 +66,7 @@ Plotting
 
 Once you have ran swift (we suggest that you use the following)
 
-    ../swift --external-gravity --stars --hydro --threads=16 keplerian_ring.yml 2>&1 | tee output.log
+    ../../../swift --external-gravity --stars --hydro --threads=16 keplerian_ring.yml 2>&1 | tee output.log
 
 there will be around 350 ```.hdf5``` files in your directory. To check out
 the results of the example use the plotting script:

examples/HydroTests/KeplerianRing/makeIC.py

@@ -35,9 +35,10 @@ class Particles(object):
     set their keplerian velocities based on their positions. These properties
     are set using the 'generationmethod' functions below.
     """
+
     def __init__(self, meta):
         self.gravitymass = meta["gravitymass"]
-        self.nparts = meta["nparts"]**2
+        self.nparts = meta["nparts"] ** 2
         self.particlemass = meta["particlemass"]
         self.softening = meta["softening"]
         self.boxsize = meta["boxsize"]
@@ -56,7 +57,6 @@ class Particles(object):
 
         return
 
-
     def calculate_masses(self):
         """
         Calculate the individual masses for the particles such that we have
@@ -66,7 +66,6 @@ class Particles(object):
         self.masses = self.densities * mass_factor
 
         return self.masses
-   
 
     def calculate_velocities(self, angle=0):
         """
@@ -75,26 +74,32 @@ class Particles(object):
 
         Please give angle in radians.
         """
-        force_modifier = np.sqrt(self.gravitymass / (self.radii**2 + self.softening**2)**(3/2)) * self.radii 
+        force_modifier = (
+            np.sqrt(
+                self.gravitymass / (self.radii ** 2 + self.softening ** 2) ** (3 / 2)
+            )
+            * self.radii
+        )
         try:
             v_x = np.sin(self.theta) * np.sin(self.phi) * np.cos(angle)
         except ValueError:
             # Phi are not yet set. Make them and then move on to v_x
             if angle:
                 raise ValueError("Unable to find phi.")
-            
+
             # If we have angle of inclination 0 we can set phi.
             self.phi = np.zeros_like(self.theta) + np.pi / 2
             v_x = np.sin(self.theta) * np.sin(self.phi) * np.cos(angle)
 
-        v_y = np.cos(self.phi) * np.sin(angle) - np.cos(self.theta) * np.sin(self.phi) * np.cos(angle)
+        v_y = np.cos(self.phi) * np.sin(angle) - np.cos(self.theta) * np.sin(
+            self.phi
+        ) * np.cos(angle)
         v_z = np.sin(self.theta) * np.sin(self.phi) * np.sin(angle)
 
         self.velocities = (force_modifier * np.array([v_x, v_y, v_z])).T
 
         return self.velocities
 
-
     def generate_ids(self):
         """
         Generate consecutive IDs from 0 based on the number of particles
@@ -104,7 +109,6 @@ class Particles(object):
 
         return self.ids
 
-
     def convert_polar_to_cartesian(self, centre_of_ring=(5, 5), boxsize=None):
         """
         Converts self.radii, self.theta to self.positions.
@@ -114,7 +118,7 @@ class Particles(object):
         if len(self.phi) == 0:
             x = self.radii * np.cos(self.theta) + centre_of_ring[0]
             y = self.radii * np.sin(self.theta) + centre_of_ring[1]
-            z = np.zeros_like(x) + boxsize/2
+            z = np.zeros_like(x) + boxsize / 2
         else:
             x = self.radii * np.cos(self.theta) * np.sin(self.phi) + centre_of_ring[0]
             y = self.radii * np.sin(self.theta) * np.sin(self.phi) + centre_of_ring[1]
@@ -122,13 +126,12 @@ class Particles(object):
                 z = self.radii * np.cos(self.phi) + centre_of_ring[2]
             except AttributeError:
                 # Just set the central z value to the middle of the box
-                z = self.radii * np.cos(self.phi) + boxsize/2
+                z = self.radii * np.cos(self.phi) + boxsize / 2
 
         self.positions = np.array([x, y, z]).T
 
         return self.positions
 
-
     def convert_cartesian_to_polar(self, centre_of_ring=(5, 5)):
         """
         Converts self.positions to self.radii, self.theta, and self.phi.
@@ -142,9 +145,9 @@ class Particles(object):
         except AttributeError:
             raise AttributeError("Invalid array for centre_of_ring provided.")
 
-        xsquare = x*x
-        ysquare = y*y
-        zsquare = z*z
+        xsquare = x * x
+        ysquare = y * y
+        zsquare = z * z
 
         r = np.sqrt(xsquare + ysquare + zsquare)
 
@@ -157,8 +160,7 @@ class Particles(object):
         theta = theta_for_unmasked + mask * 0
 
         # Thankfully we have already ensured that r != 0.
-        phi = np.arccos(z/r)
-
+        phi = np.arccos(z / r)
 
         self.radii = r
         self.theta = theta
@@ -166,7 +168,6 @@ class Particles(object):
 
         return r, theta, phi
 
-
     def wiggle_positions(self, tol=1e-3):
         """
         'Wiggle' the positions to avoid precision issues.
@@ -177,32 +178,28 @@ class Particles(object):
 
         return
 
-
     def exclude_particles(self, range):
         """
         Exclude all particles that are *not* within range (of radii).
         """
-        mask = np.logical_or(
-            self.radii < range[0],
-            self.radii > range[1],
-        )
+        mask = np.logical_or(self.radii < range[0], self.radii > range[1])
 
-        x  = np.ma.array(self.positions[:, 0], mask=mask).compressed()
-        y  = np.ma.array(self.positions[:, 1], mask=mask).compressed()
-        z  = np.ma.array(self.positions[:, 2], mask=mask).compressed()
+        x = np.ma.array(self.positions[:, 0], mask=mask).compressed()
+        y = np.ma.array(self.positions[:, 1], mask=mask).compressed()
+        z = np.ma.array(self.positions[:, 2], mask=mask).compressed()
 
         self.positions = np.array([x, y, z]).T
-        
+
         try:
-            v_x  = np.ma.array(self.velocities[:, 0], mask=mask).compressed()
-            v_y  = np.ma.array(self.velocities[:, 1], mask=mask).compressed()
-            v_z  = np.ma.array(self.velocities[:, 2], mask=mask).compressed()
+            v_x = np.ma.array(self.velocities[:, 0], mask=mask).compressed()
+            v_y = np.ma.array(self.velocities[:, 1], mask=mask).compressed()
+            v_z = np.ma.array(self.velocities[:, 2], mask=mask).compressed()
 
             self.velocities = np.array([v_x, v_y, v_z])
         except IndexError:
             # We haven't filled them yet anyway...
             pass
-    
+
         try:
             self.ids = np.ma.array(self.ids, mask=mask).compressed()
         except np.ma.core.MaskError:
@@ -218,8 +215,8 @@ class Particles(object):
         # QSP Fix has modified us, so first we need to chop off extras.
         # Then, as they are all the same, we don't need to remove 'specific'
         # values, and so we can just chop off the ends.
-        self.smoothing = self.smoothing[:len(x)]
-        self.internalenergy = self.internalenergy[:len(x)]
+        self.smoothing = self.smoothing[: len(x)]
+        self.internalenergy = self.internalenergy[: len(x)]
 
         try:
             self.masses = np.ma.array(self.masses, mask=mask).compressed()
@@ -240,7 +237,6 @@ class Particles(object):
 
         return
 
-
     def tilt_particles(self, angle, center=(5, 5, 5)):
         """
         Tilts the particles around the x-axis around the point given.
@@ -250,36 +246,23 @@ class Particles(object):
         angle_radians = angle * np.pi / 180
         rotation_matrix = np.array(
             [
-                [  np.cos(angle_radians), 0, np.sin(angle_radians)  ],
-                [            0,           1,          0             ],
-                [ -np.sin(angle_radians), 0, np.cos(angle_radians)  ]
+                [np.cos(angle_radians), 0, np.sin(angle_radians)],
+                [0, 1, 0],
+                [-np.sin(angle_radians), 0, np.cos(angle_radians)],
             ]
         )
 
         self.positions = (
-            (
-                np.matmul(
-                    rotation_matrix,
-                    (self.positions - np.array(center)).T
-                )
-            ).T
+            (np.matmul(rotation_matrix, (self.positions - np.array(center)).T)).T
         ) + np.array(center)
 
-        self.velocities = (
-            (
-                np.matmul(
-                    rotation_matrix,
-                    (self.velocities).T
-                )
-            ).T
-        )
-        
+        self.velocities = (np.matmul(rotation_matrix, (self.velocities).T)).T
+
         self.convert_cartesian_to_polar(center)
 
         return
 
-
-    def save_to_gadget(self, filename, boxsize=10.):
+    def save_to_gadget(self, filename, boxsize=10.0):
         """
         Save the particle data to a GADGET .hdf5 file.
 
@@ -293,23 +276,15 @@ class Particles(object):
                 np_total=np.array([self.nparts, 0, 0, 0, 0, 0]),
                 np_total_hw=np.array([0, 0, 0, 0, 0, 0]),
                 other={
-                    "MassTable" : np.array([self.particlemass, 0, 0, 0, 0, 0]),
-                    "Time" : 0,
-                }
+                    "MassTable": np.array([self.particlemass, 0, 0, 0, 0, 0]),
+                    "Time": 0,
+                },
             )
 
-            wg.write_runtime_pars(
-                handle,
-                periodic_boundary=1,
-            )
+            wg.write_runtime_pars(handle, periodic_boundary=1)
 
             wg.write_units(
-                handle,
-                current=1.,
-                length=1.,
-                mass=1,
-                temperature=1.,
-                time=1.,
+                handle, current=1.0, length=1.0, mass=1, temperature=1.0, time=1.0
             )
 
             wg.write_block(
@@ -332,11 +307,11 @@ def __sigma(r):
     Density distribution of the ring, this comes directly from Hopkins 2015.
     """
     if r < 0.5:
-        return 0.01 + (r/0.5)**3
+        return 0.01 + (r / 0.5) ** 3
     elif r <= 2:
         return 1.01
     else:
-        return 0.01 + (1 + (r-2)/0.1)**(-3)
+        return 0.01 + (1 + (r - 2) / 0.1) ** (-3)
 
 
 # This is required because of the if, else statement above that does not
@@ -344,7 +319,7 @@ def __sigma(r):
 sigma = np.vectorize(__sigma)
 
 
-def generate_theta(r, theta_initial=0.):
+def generate_theta(r, theta_initial=0.0):
     """
     Generate the theta associated with the particles based on their radii.
     This uses the method from The effect of Poisson noise on SPH calculations,
@@ -370,7 +345,7 @@ def generate_theta(r, theta_initial=0.):
     # Need to do this awful loop because each entry relies on the one before it.
     # Unless someone knows how to do this in numpy.
     for i in range(len(d_theta_i)):  # first is theta_initial
-        theta_i[i+1] = theta_i[i] + d_theta_i[i]
+        theta_i[i + 1] = theta_i[i] + d_theta_i[i]
 
     return theta_i
 
@@ -425,7 +400,7 @@ def QSP_fix(r_i, theta_i):
     # We want to have all particles in each circle a fixed radius away from the
     # centre, as well as having even spacing between each particle. The final
     # ring may be a bit dodgy still, but we will see.
-    
+
     r_i_fixed = []
     theta_i_fixed = []
 
@@ -435,7 +410,7 @@ def QSP_fix(r_i, theta_i):
         theta_initial = circle[1][0]
 
         theta_sep = 2 * np.pi / n_particles
-        
+
         theta = [t * theta_sep for t in range(n_particles)]
         radii = [radius] * n_particles
 
@@ -451,24 +426,26 @@ def gen_particles_grid(meta):
     """
     particles = Particles(meta)
     range = (0, meta["boxsize"])
-    centre_of_ring = [meta["boxsize"]/2.] * 3
+    centre_of_ring = [meta["boxsize"] / 2.0] * 3
 
     # Because we are using a uniform grid we actually use the same x and y
     # range for the initial particle setup.
-    step = (range[1] - range[0])/meta["nparts"]
-    
+    step = (range[1] - range[0]) / meta["nparts"]
+
     x_values = np.arange(0, range[1] - range[0], step)
 
     # These are 2d arrays which isn't actually that helpful.
     x, y = np.meshgrid(x_values, x_values)
-    x = x.flatten() + centre_of_ring[0] - (range[1] - range[0])/2
-    y = y.flatten() + centre_of_ring[1] - (range[1] - range[0])/2
-    z = np.zeros_like(x) + meta["boxsize"]/2
+    x = x.flatten() + centre_of_ring[0] - (range[1] - range[0]) / 2
+    y = y.flatten() + centre_of_ring[1] - (range[1] - range[0]) / 2
+    z = np.zeros_like(x) + meta["boxsize"] / 2
 
     particles.positions = np.array([x, y, z]).T
-    particles.radii = np.sqrt((x - centre_of_ring[0])**2 + (y - centre_of_ring[1])**2)
+    particles.radii = np.sqrt(
+        (x - centre_of_ring[0]) ** 2 + (y - centre_of_ring[1]) ** 2
+    )
     particles.theta = np.arctan2(y - centre_of_ring[1], x - centre_of_ring[0])
-    particles.exclude_particles((particles.softening, 100.))
+    particles.exclude_particles((particles.softening, 100.0))
 
     particles.densities = sigma(particles.radii)
     particles.calculate_velocities()
@@ -488,23 +465,27 @@ def gen_particles_spiral(meta):
     object. Based on Cartwright, Stamatellos & Whitworth (2009).
     """
     particles = Particles(meta)
-    centre_of_ring = (meta["boxsize"]/2., meta["boxsize"]/2., meta["boxsize"]/2.)
-    max_r = meta["boxsize"]/2.
+    centre_of_ring = (
+        meta["boxsize"] / 2.0,
+        meta["boxsize"] / 2.0,
+        meta["boxsize"] / 2.0,
+    )
+    max_r = meta["boxsize"] / 2.0
 
-    m = (np.arange(particles.nparts) + 0.5)/particles.nparts
+    m = (np.arange(particles.nparts) + 0.5) / particles.nparts
     r = max_r * m
     theta = generate_theta(r)
 
     particles.radii, particles.theta = QSP_fix(r, theta)
     # We must do this afterwards as QSP does not always return the same number of parts.
-    phi = np.zeros_like(particles.theta) + np.pi/2
+    phi = np.zeros_like(particles.theta) + np.pi / 2
     particles.phi = phi
 
     particles.convert_polar_to_cartesian(centre_of_ring, meta["boxsize"])
     particles.nparts = len(particles.radii)
 
-    particles.exclude_particles((particles.softening, 100.))
-    
+    particles.exclude_particles((particles.softening, 100.0))
+
     # This way of doing densities does not work for different sized patches.
     # Therefore we need to weight by effecitve area.
     dtheta = (particles.theta[1:] - particles.theta[:-1]) / 2
@@ -521,7 +502,6 @@ def gen_particles_spiral(meta):
     if meta["angle"]:
         particles.tilt_particles(meta["angle"], centre_of_ring)
 
-
     return particles
 
 
@@ -533,25 +513,29 @@ def gen_particles_gaussian(meta):
     This generation function uses a Gaussian PDF.
     """
     particles = Particles(meta)
-    centre_of_ring = (meta["boxsize"]/2., meta["boxsize"]/2., meta["boxsize"]/2.)
-    max_r = meta["boxsize"]/2.
+    centre_of_ring = (
+        meta["boxsize"] / 2.0,
+        meta["boxsize"] / 2.0,
+        meta["boxsize"] / 2.0,
+    )
+    max_r = meta["boxsize"] / 2.0
 
-    m = (np.arange(particles.nparts) + 0.5)/particles.nparts
-    error_function = erfinv(2*m - 1)
+    m = (np.arange(particles.nparts) + 0.5) / particles.nparts
+    error_function = erfinv(2 * m - 1)
     r = 2 + 0.5 * error_function
 
     theta = generate_theta(r)
 
     particles.radii, particles.theta = QSP_fix(r, theta)
     # We must do this afterwards as QSP does not always return the same number of parts.
-    phi = np.zeros_like(particles.theta) + np.pi/2
+    phi = np.zeros_like(particles.theta) + np.pi / 2
     particles.phi = phi
 
     particles.convert_polar_to_cartesian(centre_of_ring, meta["boxsize"])
     particles.nparts = len(particles.radii)
 
-    particles.exclude_particles((particles.softening, 100.))
-    
+    particles.exclude_particles((particles.softening, 100.0))
+
     # This way of doing densities does not work for different sized patches.
     particles.densities = np.zeros_like(particles.radii)
     particles.calculate_velocities()
@@ -562,7 +546,6 @@ def gen_particles_gaussian(meta):
     if meta["angle"]:
         particles.tilt_particles(meta["angle"], centre_of_ring)
 
-
     return particles
 
 
@@ -585,7 +568,7 @@ if __name__ == "__main__":
              GM for the central point mass. Default: 1.
              """,
         required=False,
-        default=1.,
+        default=1.0,
     )
 
     PARSER.add_argument(
@@ -619,7 +602,7 @@ if __name__ == "__main__":
              Maximum mass of the gas particles. Default: 10.
              """,
         required=False,
-        default=10.,
+        default=10.0,
     )
 
     PARSER.add_argument(
@@ -644,7 +627,7 @@ if __name__ == "__main__":
         required=False,
         default=-1,
     )
-    
+
     PARSER.add_argument(
         "-i",
         "--internalenergy",
@@ -654,7 +637,7 @@ if __name__ == "__main__":
              SPH for details). Default: 0.015.
              """,
         required=False,
-        default=0.015
+        default=0.015,
     )
 
     PARSER.add_argument(
@@ -678,7 +661,7 @@ if __name__ == "__main__":
              The box size. Default: 10.
              """,
         required=False,
-        default=10.
+        default=10.0,
     )
 
     PARSER.add_argument(
@@ -691,10 +674,9 @@ if __name__ == "__main__":
              particles.
              """,
         required=False,
-        default=0.
+        default=0.0,
     )
 
-
     ### --- ### --- Argument Parsing --- ### --- ###
 
     ARGS = vars(PARSER.parse_args())
@@ -713,13 +695,11 @@ if __name__ == "__main__":
         )
         exit(1)
 
-
     if ARGS["smoothing"] == -1:
         smoothing = float(ARGS["boxsize"]) / int(ARGS["nparts"])
     else:
         smoothing = float(ARGS["smoothing"])
 
-
     META = {
         "gravitymass": float(ARGS["gravitymass"]),
         "nparts": int(ARGS["nparts"]),
@@ -728,13 +708,9 @@ if __name__ == "__main__":
         "softening": float(ARGS["softening"]),
         "internalenergy": float(ARGS["internalenergy"]),
         "boxsize": float(ARGS["boxsize"]),
-        "angle" : float(ARGS["angle"])
+        "angle": float(ARGS["angle"]),
     }
 
     PARTICLES = gen_particles(META)
 
-    PARTICLES.save_to_gadget(
-        filename=ARGS["filename"],
-        boxsize=ARGS["boxsize"],
-    )
-
+    PARTICLES.save_to_gadget(filename=ARGS["filename"], boxsize=ARGS["boxsize"])

examples/HydroTests/KeplerianRing/make_movie.py

@@ -28,6 +28,7 @@
 """
 
 import matplotlib
+
 matplotlib.use("Agg")
 
 
@@ -46,36 +47,34 @@ def get_colours(numbers, norm=None, cm_name="viridis"):
     cmap = matplotlib.cm.get_cmap(cm_name)
 
     return cmap(norm(numbers)), norm
-        
+
 
 def load_data(filename, silent=True, extra=None, logextra=True, exclude=None):
     if not silent:
         print(f"Loading data from {filename}")
 
     with h5.File(filename, "r") as file_handle:
-        coords = file_handle['PartType0']['Coordinates'][...]
-        time = float(file_handle['Header'].attrs['Time'])
-        boxsize = file_handle['Header'].attrs['BoxSize']
+        coords = file_handle["PartType0"]["Coordinates"][...]
+        time = float(file_handle["Header"].attrs["Time"])
+        boxsize = file_handle["Header"].attrs["BoxSize"]
 
         # For some old runs we have z=0
         if np.sum(coords[:, 2]) == 0:
-            centre_of_box = np.array(list(boxsize[:2]/2) + [0])
+            centre_of_box = np.array(list(boxsize[:2] / 2) + [0])
         else:
-            centre_of_box = boxsize/2
+            centre_of_box = boxsize / 2
 
-        
         if exclude is not None:
             distance_from_centre = coords - centre_of_box
             r2 = np.sum(distance_from_centre * distance_from_centre, 1)
             mask = r2 < (exclude * exclude)
 
-            masked = np.ma.array(coords, mask=np.array([mask]*3))
+            masked = np.ma.array(coords, mask=np.array([mask] * 3))
 
             coords = masked.compressed()
-        
 
         if extra is not None:
-            other = file_handle['PartType0'][extra][...]
+            other = file_handle["PartType0"][extra][...]
             if exclude is not None:
                 masked = np.ma.array(other, mask=mask)
 
@@ -90,7 +89,7 @@ def load_data(filename, silent=True, extra=None, logextra=True, exclude=None):
 
 
 def rms(x):
-    return np.sqrt(sum(x**2))
+    return np.sqrt(sum(x ** 2))
 
 
 def rotation_velocity_at_r(r, params):
@@ -102,7 +101,7 @@ def rotation_velocity_at_r(r, params):
 
     unit_length = float(params[r"InternalUnitSystem:UnitLength_in_cgs"])
 
-    if unit_length != 1.:
+    if unit_length != 1.0:
         print(f"Your unit length: {unit_length}")
         raise InternalUnitSystemError(
             "This function is only created to handle CGS units."
@@ -111,7 +110,7 @@ def rotation_velocity_at_r(r, params):
     central_mass = float(params["PointMassPotential:mass"])
     G = 6.67408e-8
 
-    v = np.sqrt( G * central_mass / r)
+    v = np.sqrt(G * central_mass / r)
 
     return v
 
@@ -123,25 +122,25 @@ def get_rotation_period_at_r(r, params):
     """
     v = rotation_velocity_at_r(r, params)
 
-    return 2*np.pi / v
+    return 2 * np.pi / v
 
 
 def get_metadata(filename, r=1):
     """ The metadata should be extracted from the first snapshot. """
     with h5.File(filename, "r") as file_handle:
-        header = file_handle['Header'].attrs
-        code = file_handle['Code'].attrs
-        hydro = file_handle['HydroScheme'].attrs
-        params = file_handle['Parameters'].attrs
+        header = file_handle["Header"].attrs
+        code = file_handle["Code"].attrs
+        hydro = file_handle["HydroScheme"].attrs
+        params = file_handle["Parameters"].attrs
 
         period = get_rotation_period_at_r(r, params)
 
         return_values = {
-            "header" : dict(header),
-            "code" : dict(code),
-            "period" : float(period),
-            "hydro" : dict(hydro),
-            "params" : dict(params)
+            "header": dict(header),
+            "code": dict(code),
+            "period": float(period),
+            "hydro": dict(hydro),
+            "params": dict(params),
         }
 
     return return_values
@@ -155,12 +154,8 @@ def plot_single(number, scatter, text, metadata, ax, extra=None, norm=None):
     else:
         coordinates, time = load_data(filename)
 
-
     text.set_text(
-        "Time: {:1.2f} | Rotations {:1.2f}".format(
-            time,
-            time/metadata['period'],
-        )
+        "Time: {:1.2f} | Rotations {:1.2f}".format(time, time / metadata["period"])
     )
 
     data = coordinates[:, 0:2]
@@ -170,7 +165,7 @@ def plot_single(number, scatter, text, metadata, ax, extra=None, norm=None):
         colours, _ = get_colours(other, norm)
         scatter.set_color(colours)
 
-    return scatter,
+    return (scatter,)
 
 
 if __name__ == "__main__":
@@ -193,8 +188,8 @@ if __name__ == "__main__":
              still plotted -- they are just excluded for the purposes
              of colourmapping. Default: 1 simulation unit.
              """,
-        default=1.,
-        required=False
+        default=1.0,
+        required=False,
     )
 
     parser.add_argument(
@@ -206,7 +201,7 @@ if __name__ == "__main__":
              Default: don't use a colourmap. (Much faster).
              """,
         required=False,
-        default=None
+        default=None,
     )
 
     parser.add_argument(
@@ -219,7 +214,7 @@ if __name__ == "__main__":
         required=False,
         default=2.5,
     )
-    
+
     args = vars(parser.parse_args())
 
     # Look for the number of files in the directory.
@@ -233,12 +228,18 @@ if __name__ == "__main__":
         if i > 10000:
             break
 
-
-   
     # Now deal with the colourmapping (if present)
     if args["cmap"] is not None:
-        _, _, numbers0 = load_data("keplerian_ring_0000.hdf5", extra=args["cmap"], exclude=float(args["exclude_central"]))
-        _, _, numbersend = load_data("keplerian_ring_{:04d}.hdf5".format(i-1), extra=args["cmap"], exclude=float(args["exclude_central"]))
+        _, _, numbers0 = load_data(
+            "keplerian_ring_0000.hdf5",
+            extra=args["cmap"],
+            exclude=float(args["exclude_central"]),
+        )
+        _, _, numbersend = load_data(
+            "keplerian_ring_{:04d}.hdf5".format(i - 1),
+            extra=args["cmap"],
+            exclude=float(args["exclude_central"]),
+        )
         vmax = max([numbers0.max(), numbersend.max()])
         vmin = min([numbers0.min(), numbersend.min()])
 
@@ -246,19 +247,18 @@ if __name__ == "__main__":
     else:
         norm = None
 
-
     # Initial plot setup
 
     metadata = get_metadata("keplerian_ring_0000.hdf5")
-    n_particle = metadata['header']['NumPart_Total'][0]
+    n_particle = metadata["header"]["NumPart_Total"][0]
 
     fig = plt.figure(figsize=(8, 8))
     ax = fig.add_subplot(111)
 
-    scatter = ax.scatter([0]*n_particle, [0]*n_particle, s=0.5, marker="o")
+    scatter = ax.scatter([0] * n_particle, [0] * n_particle, s=0.5, marker="o")
     diff = float(args["max"])
-    left = metadata['header']['BoxSize'][0]/2 - diff
-    right = metadata['header']['BoxSize'][0]/2 + diff
+    left = metadata["header"]["BoxSize"][0] / 2 - diff
+    right = metadata["header"]["BoxSize"][0] / 2 + diff
     ax.set_xlim(left, right)
     ax.set_ylim(left, right)
 
@@ -266,45 +266,34 @@ if __name__ == "__main__":
     time_text = ax.text(
         offset + left,
         offset + left,
-        "Time: {:1.2f} | Rotations {:1.2f}".format(
-            0,
-            0/metadata['period'],
-        )
+        "Time: {:1.2f} | Rotations {:1.2f}".format(0, 0 / metadata["period"]),
     )
 
     ax.text(
         offset + left,
-        right-offset-0.35,
+        right - offset - 0.35,
         "Code: {} {} | {} {} \nHydro {}\n$\eta$={:1.4f}".format(
-            metadata['code']['Git Branch'].decode("utf-8"),
-            metadata['code']['Git Revision'].decode("utf-8"),
-            metadata['code']['Compiler Name'].decode("utf-8"),
-            metadata['code']['Compiler Version'].decode("utf-8"),
-            metadata['hydro']['Scheme'].decode("utf-8"),
-            metadata['hydro']['Kernel eta'][0],
-        )
+            metadata["code"]["Git Branch"].decode("utf-8"),
+            metadata["code"]["Git Revision"].decode("utf-8"),
+            metadata["code"]["Compiler Name"].decode("utf-8"),
+            metadata["code"]["Compiler Version"].decode("utf-8"),
+            metadata["hydro"]["Scheme"].decode("utf-8"),
+            metadata["hydro"]["Kernel eta"][0],
+        ),
     )
 
     ax.set_title("Keplerian Ring Test")
     ax.set_xlabel("$x$ position")
     ax.set_ylabel("$y$ position")
 
-    
     anim = anim.FuncAnimation(
         fig,
         plot_single,
         tqdm(np.arange(i)),
-        fargs = [
-            scatter,
-            time_text,
-            metadata,
-            ax,
-            args["cmap"],
-            norm
-        ],
+        fargs=[scatter, time_text, metadata, ax, args["cmap"], norm],
         interval=50,
         repeat_delay=3000,
         blit=True,
     )
 
-    anim.save("keplerian_ring.mp4", dpi=int(640/8))
+    anim.save("keplerian_ring.mp4", dpi=int(640 / 8))