diff --git a/examples/Planetary/JupiterLikePlanet/make_init_cond.py b/examples/Planetary/JupiterLikePlanet/make_init_cond.py
index af6cc518cdaac2a72e17f7ca7ecf774ed4bd05d6..e32513fcb1d946e878934f6628968080eb618afe 100755
--- a/examples/Planetary/JupiterLikePlanet/make_init_cond.py
+++ b/examples/Planetary/JupiterLikePlanet/make_init_cond.py
@@ -1,7 +1,7 @@
 ################################################################################
 # This file is part of SWIFT.
 # Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
-#		        2024 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+# 		        2024 Jacob Kegerreis (jacob.kegerreis@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
@@ -24,7 +24,7 @@ import h5py
 import woma
 
 # Number of particles
-N = 10 ** 7
+N = 10**7
 N_label = "n%d" % (10 * np.log10(N))
 
 # Earth units
@@ -41,16 +41,14 @@ target_prof = woma.Planet(
     A1_M_layer=[10 * M_E, 298 * M_E],
     P_s=1e5,
     T_s=165,
-    num_prof= 10000,
+    num_prof=10000,
 )
 
 # Load material tables
-woma.load_eos_tables(
-    np.unique(target_prof.A1_mat_layer)
-)
+woma.load_eos_tables(np.unique(target_prof.A1_mat_layer))
 
 # Compute profiles
-target_prof.gen_prof_L2_find_R_R1_given_M1_M2(R_min = 10.7*R_E, R_max = 11.2*R_E)
+target_prof.gen_prof_L2_find_R_R1_given_M1_M2(R_min=10.7 * R_E, R_max=11.2 * R_E)
 
 # Save profile data
 target_prof.save("demo_target_profile.hdf5")
diff --git a/examples/Planetary/JupiterLikePlanet/plot_profiles.py b/examples/Planetary/JupiterLikePlanet/plot_profiles.py
index da43aa9df034034bbb4781ae6953e004da317357..8bba4541fd0ac81bf831adde6efbdeaa3bd7b19a 100755
--- a/examples/Planetary/JupiterLikePlanet/plot_profiles.py
+++ b/examples/Planetary/JupiterLikePlanet/plot_profiles.py
@@ -1,7 +1,7 @@
 ###############################################################################
 # This file is part of SWIFT.
 # Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
-#	            2023 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+# 	            2023 Jacob Kegerreis (jacob.kegerreis@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
@@ -35,7 +35,7 @@ import h5py
 import woma
 
 # Number of particles
-N = 10 ** 7
+N = 10**7
 N_label = "n%d" % (10 * np.log10(N))
 
 # Plotting options
@@ -62,7 +62,7 @@ def plot_profile_and_particles(profile, A1_r, A1_rho):
     ax.plot(profile.A1_r / R_E, profile.A1_rho)
 
     # Particles
-    ax.scatter(A1_r / R_E, A1_rho, c="k", marker=".", s=1 ** 2)
+    ax.scatter(A1_r / R_E, A1_rho, c="k", marker=".", s=1**2)
 
     ax.set_xlim(0, None)
     ax.set_xlabel(r"Radial distance ($R_\oplus$)")
@@ -93,7 +93,7 @@ if __name__ == "__main__":
                 np.array(f["PartType0/Coordinates"][()])
                 - 0.5 * f["Header"].attrs["BoxSize"]
             ) * file_to_SI.l
-            A1_r = np.sqrt(np.sum(A2_pos ** 2, axis=1))
+            A1_r = np.sqrt(np.sum(A2_pos**2, axis=1))
             A1_rho = np.array(f["PartType0/Densities"][()]) * file_to_SI.rho
 
         # Plot the data
diff --git a/examples/Planetary/JupiterLikePlanet/plot_snapshots.py b/examples/Planetary/JupiterLikePlanet/plot_snapshots.py
index 2bcada8746939f1dd16432c6354c8013904e6fb8..a6db7c12dd22c5d9544ac1f723662a9608d576b6 100755
--- a/examples/Planetary/JupiterLikePlanet/plot_snapshots.py
+++ b/examples/Planetary/JupiterLikePlanet/plot_snapshots.py
@@ -1,7 +1,7 @@
 ###############################################################################
 # This file is part of SWIFT.
 # Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
-#	            2023 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+# 	            2023 Jacob Kegerreis (jacob.kegerreis@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
@@ -28,7 +28,7 @@ import h5py
 import woma
 
 # Number of particles
-N = 10 ** 7
+N = 10**7
 N_label = "n%d" % (10 * np.log10(N))
 
 # Plotting options
@@ -47,7 +47,7 @@ Di_mat_colour = {"AQUA": "orangered", "CD21_HHe": "gold"}
 Di_id_colour = {woma.Di_mat_id[mat]: colour for mat, colour in Di_mat_colour.items()}
 
 # Scale point size with resolution
-size = (1 * np.cbrt(10 ** 6 / N)) ** 2
+size = (1 * np.cbrt(10**6 / N)) ** 2
 
 
 def load_snapshot(filename):
diff --git a/examples/Planetary/KelvinHelmholtz_EarthLike_3D/makeIC.py b/examples/Planetary/KelvinHelmholtz_EarthLike_3D/makeIC.py
index 4f3fda91833fb45ffc8bc275235ab9e64fb615ed..7745cff97a5280bfaa6cce87ab9cd3464d61846a 100644
--- a/examples/Planetary/KelvinHelmholtz_EarthLike_3D/makeIC.py
+++ b/examples/Planetary/KelvinHelmholtz_EarthLike_3D/makeIC.py
@@ -1,22 +1,22 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
- #               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) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+#               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/>.
+#
+##############################################################################
 
 import h5py
 import numpy as np
@@ -24,26 +24,26 @@ import numpy as np
 # Generates a swift IC file for the Kelvin-Helmholtz test in a periodic box
 
 # Constants
-R_earth = 6371000    # Earth radius
+R_earth = 6371000  # Earth radius
 
 # Parameters
-N2_l     = 128       # Particles along one edge in the low-density region
-N2_depth = 18        # Particles in z direction in low-density region
-matID1 = 402         # Central region material ID: ANEOS Fe85Si15
-matID2 = 400         # Outskirts material ID: ANEOS forsterite
-P1    = 1.2e+11      # Central region pressure
-P2    = 1.2e+11      # Outskirts pressure
-u1    = 4069874      # Central region specific internal energy
-u2    = 9899952      # Outskirts specific internal energy
-rho1_approx =  10000 # Central region density. Readjusted later
-rho2 =  5000         # Outskirts density
-boxsize_l = R_earth        # size of simulation box in x and y dimension
-v1   = boxsize_l / 10000  # Central region velocity
-v2   = -boxsize_l / 10000 # Outskirts velocity
-boxsize_depth = boxsize_l * N2_depth / N2_l # size of simulation box in z dimension
+N2_l = 128  # Particles along one edge in the low-density region
+N2_depth = 18  # Particles in z direction in low-density region
+matID1 = 402  # Central region material ID: ANEOS Fe85Si15
+matID2 = 400  # Outskirts material ID: ANEOS forsterite
+P1 = 1.2e11  # Central region pressure
+P2 = 1.2e11  # Outskirts pressure
+u1 = 4069874  # Central region specific internal energy
+u2 = 9899952  # Outskirts specific internal energy
+rho1_approx = 10000  # Central region density. Readjusted later
+rho2 = 5000  # Outskirts density
+boxsize_l = R_earth  # size of simulation box in x and y dimension
+v1 = boxsize_l / 10000  # Central region velocity
+v2 = -boxsize_l / 10000  # Outskirts velocity
+boxsize_depth = boxsize_l * N2_depth / N2_l  # size of simulation box in z dimension
 mass = rho2 * (boxsize_l * boxsize_l * boxsize_depth) / (N2_l * N2_l * N2_depth)
 fileOutputName = "kelvin_helmholtz.hdf5"
-#---------------------------------------------------
+# ---------------------------------------------------
 
 # Start by calculating N1_l and rho1
 numPart2 = N2_l * N2_l * N2_depth
@@ -86,23 +86,27 @@ for i in range(N1_depth):
     for j in range(N1_l):
         for k in range(N1_l):
             index = i * N1_l**2 + j * N1_l + k
-            A2_coords1[index, 0] = (j / float(N1_l) + 1. / (2. * N1_l)) * boxsize_l
-            A2_coords1[index, 1] = (k / float(N1_l) + 1. / (2. * N1_l)) * boxsize_l
-            A2_coords1[index, 2] = (i / float(N1_depth) + 1. / (2. * N1_depth)) * boxsize_depth
+            A2_coords1[index, 0] = (j / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_l
+            A2_coords1[index, 1] = (k / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_l
+            A2_coords1[index, 2] = (
+                i / float(N1_depth) + 1.0 / (2.0 * N1_depth)
+            ) * boxsize_depth
 
 # Particles in the outskirts
 for i in range(N2_depth):
     for j in range(N2_l):
         for k in range(N2_l):
             index = i * N2_l**2 + j * N2_l + k
-            A2_coords2[index, 0] = (j / float(N2_l) + 1. / (2. * N2_l)) * boxsize_l
-            A2_coords2[index, 1] = (k / float(N2_l) + 1. / (2. * N2_l)) * boxsize_l
-            A2_coords2[index, 2] = (i / float(N2_depth) + 1. / (2. * N2_depth)) * boxsize_depth
+            A2_coords2[index, 0] = (j / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_l
+            A2_coords2[index, 1] = (k / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_l
+            A2_coords2[index, 2] = (
+                i / float(N2_depth) + 1.0 / (2.0 * N2_depth)
+            ) * boxsize_depth
 
 
 # Masks for the particles to be selected for the outer and inner regions
-mask1 = abs(A2_coords1[:,1] - 0.5 * boxsize_l) < 0.25 * boxsize_l
-mask2 = abs(A2_coords2[:,1] - 0.5 * boxsize_l) > 0.25 * boxsize_l
+mask1 = abs(A2_coords1[:, 1] - 0.5 * boxsize_l) < 0.25 * boxsize_l
+mask2 = abs(A2_coords2[:, 1] - 0.5 * boxsize_l) > 0.25 * boxsize_l
 
 # The positions of the particles are now selected
 # and the placement of the lattices are adjusted to give appropriate interfaces
@@ -115,16 +119,22 @@ pcl_separation_2 = np.cbrt(mass / rho2)
 boundary_separation = 0.5 * (pcl_separation_1 + pcl_separation_2)
 
 # Shift all the "inside" particles to get boundary_separation across the bottom interface
-min_y_inside = np.min(A2_coords_inside[:,1])
-max_y_outside_bot = np.max(A2_coords_outside[A2_coords_outside[:,1]-0.5 * boxsize_l < -0.25 * boxsize_l, 1])
+min_y_inside = np.min(A2_coords_inside[:, 1])
+max_y_outside_bot = np.max(
+    A2_coords_outside[A2_coords_outside[:, 1] - 0.5 * boxsize_l < -0.25 * boxsize_l, 1]
+)
 shift_distance_bot = boundary_separation - (min_y_inside - max_y_outside_bot)
-A2_coords_inside[:,1] += shift_distance_bot
+A2_coords_inside[:, 1] += shift_distance_bot
 
 # Shift the top section of the "outside" particles to get boundary_separation across the top interface
-max_y_inside = np.max(A2_coords_inside[:,1])
-min_y_outside_top = np.min(A2_coords_outside[A2_coords_outside[:,1]-0.5 * boxsize_l > 0.25 * boxsize_l, 1])
+max_y_inside = np.max(A2_coords_inside[:, 1])
+min_y_outside_top = np.min(
+    A2_coords_outside[A2_coords_outside[:, 1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1]
+)
 shift_distance_top = boundary_separation - (min_y_outside_top - max_y_inside)
-A2_coords_outside[A2_coords_outside[:,1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1] += shift_distance_top
+A2_coords_outside[
+    A2_coords_outside[:, 1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1
+] += shift_distance_top
 
 # Adjust box size in y direction based on the shifting of the lattices.
 new_box_y = boxsize_l + shift_distance_top
@@ -142,14 +152,16 @@ A1_ids = np.linspace(1, numPart, numPart)
 
 # Finally add the velocity perturbation
 vel_perturb_factor = 0.01 * (v1 - v2)
-A2_vel[:,1] = vel_perturb_factor * np.sin(2 * np.pi * A2_coords[:,0] / (0.5 * boxsize_l))
+A2_vel[:, 1] = vel_perturb_factor * np.sin(
+    2 * np.pi * A2_coords[:, 0] / (0.5 * boxsize_l)
+)
 
 # Write ICs to file
 with h5py.File(fileOutputName, "w") as f:
     # Header
     grp = f.create_group("/Header")
     grp.attrs["BoxSize"] = [boxsize_l, new_box_y, boxsize_depth]
-    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
@@ -158,29 +170,29 @@ with h5py.File(fileOutputName, "w") as f:
     grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
     grp.attrs["Dimension"] = 3
 
-    #Units
+    # Units
     grp = f.create_group("/Units")
-    grp.attrs["Unit length in cgs (U_L)"] = 100.
-    grp.attrs["Unit mass in cgs (U_M)"] = 1000.
-    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)"] = 100.0
+    grp.attrs["Unit mass in cgs (U_M)"] = 1000.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 = f.create_group("/PartType0")
-    ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+    ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
     ds[()] = A2_coords
-    ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+    ds = grp.create_dataset("Velocities", (numPart, 3), "f")
     ds[()] = A2_vel
-    ds = grp.create_dataset('Masses', (numPart, 1), 'f')
-    ds[()] = A1_m.reshape((numPart,1))
-    ds = grp.create_dataset('Density', (numPart,1), 'f')
-    ds[()] = A1_rho.reshape((numPart,1))
-    ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-    ds[()] = A1_h.reshape((numPart,1))
-    ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-    ds[()] = A1_u.reshape((numPart,1))
-    ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
-    ds[()] = A1_ids.reshape((numPart,1))
-    ds = grp.create_dataset('MaterialIDs', (numPart,1), 'i')
-    ds[()] = A1_mat.reshape((numPart,1))
+    ds = grp.create_dataset("Masses", (numPart, 1), "f")
+    ds[()] = A1_m.reshape((numPart, 1))
+    ds = grp.create_dataset("Density", (numPart, 1), "f")
+    ds[()] = A1_rho.reshape((numPart, 1))
+    ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+    ds[()] = A1_h.reshape((numPart, 1))
+    ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+    ds[()] = A1_u.reshape((numPart, 1))
+    ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+    ds[()] = A1_ids.reshape((numPart, 1))
+    ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
+    ds[()] = A1_mat.reshape((numPart, 1))
diff --git a/examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotModeGrowth.py b/examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotModeGrowth.py
index 5ff9211bd4031c14463c75d7ce237f6b304ee5bd..3785996dc39a1f31d92e1ed9351b5ff20b06f875 100644
--- a/examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotModeGrowth.py
+++ b/examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotModeGrowth.py
@@ -28,6 +28,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def calculate_mode_growth(snap, vy_init_amp, wavelength):
 
     # Load data
@@ -40,46 +41,72 @@ def calculate_mode_growth(snap, vy_init_amp, wavelength):
         A1_h = f["/PartType0/SmoothingLengths"][:] / boxsize_l
 
     # Adjust positions
-    A2_pos[:,0] += 0.5
-    A2_pos[:,1] += 0.5
+    A2_pos[:, 0] += 0.5
+    A2_pos[:, 1] += 0.5
 
     # Masks to select the upper and lower halfs of the simulations
-    mask_up = A2_pos[:,1] >= 0.5
-    mask_down = A2_pos[:,1] < 0.5
+    mask_up = A2_pos[:, 1] >= 0.5
+    mask_down = A2_pos[:, 1] < 0.5
 
     # McNally et al. 2012 Eqn. 10
     s = np.empty(len(A1_h))
-    s[mask_down] = A2_vel[mask_down, 1] * A1_h[mask_down]**3 * np.sin(2 * np.pi * A2_pos[mask_down, 0] / wavelength) * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    s[mask_up] = A2_vel[mask_up, 1] * A1_h[mask_up]**3 * np.sin(2 * np.pi * A2_pos[mask_up, 0] / wavelength) * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    s[mask_down] = (
+        A2_vel[mask_down, 1]
+        * A1_h[mask_down] ** 3
+        * np.sin(2 * np.pi * A2_pos[mask_down, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
+    )
+    s[mask_up] = (
+        A2_vel[mask_up, 1]
+        * A1_h[mask_up] ** 3
+        * np.sin(2 * np.pi * A2_pos[mask_up, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    )
 
     # McNally et al. 2012 Eqn. 11
     c = np.empty(len(A1_h))
-    c[mask_down] = A2_vel[mask_down, 1] * A1_h[mask_down]**3 * np.cos(2 * np.pi * A2_pos[mask_down, 0] / wavelength) * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    c[mask_up] = A2_vel[mask_up, 1] * A1_h[mask_up]**3 * np.cos(2 * np.pi * A2_pos[mask_up, 0] / wavelength) * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    c[mask_down] = (
+        A2_vel[mask_down, 1]
+        * A1_h[mask_down] ** 3
+        * np.cos(2 * np.pi * A2_pos[mask_down, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
+    )
+    c[mask_up] = (
+        A2_vel[mask_up, 1]
+        * A1_h[mask_up] ** 3
+        * np.cos(2 * np.pi * A2_pos[mask_up, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    )
 
     # McNally et al. 2012 Eqn. 12
     d = np.empty(len(A1_h))
-    d[mask_down] = A1_h[mask_down]**3 * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    d[mask_up] = A1_h[mask_up]**3 * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    d[mask_down] = A1_h[mask_down] ** 3 * np.exp(
+        -2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength
+    )
+    d[mask_up] = A1_h[mask_up] ** 3 * np.exp(
+        -2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength
+    )
 
     # McNally et al. 2012 Eqn. 13
-    M = 2 * np.sqrt((np.sum(s)/ np.sum(d))**2 + (np.sum(c)/ np.sum(d))**2)
+    M = 2 * np.sqrt((np.sum(s) / np.sum(d)) ** 2 + (np.sum(c) / np.sum(d)) ** 2)
 
     return M
 
 
-
 if __name__ == "__main__":
 
     # Simulation paramerters for nomralisation of mode growth
-    vy_init_amp = 2e-6 # Initial amplitude of y velocity perturbation in units of the boxsize per second
-    wavelength = 0.5   # wavelength of initial perturbation in units of the boxsize
+    vy_init_amp = 2e-6  # Initial amplitude of y velocity perturbation in units of the boxsize per second
+    wavelength = 0.5  # wavelength of initial perturbation in units of the boxsize
 
     # Use Gridspec to set up figure
     n_gs_ax = 40
     ax_len = 5
     fig = plt.figure(figsize=(9, 6))
-    gs = mpl.gridspec.GridSpec(nrows=40, ncols=60,)
+    gs = mpl.gridspec.GridSpec(
+        nrows=40,
+        ncols=60,
+    )
     ax = plt.subplot(gs[:, :])
 
     # Snapshots and corresponding times
@@ -98,8 +125,8 @@ if __name__ == "__main__":
     ax.set_ylabel(r"$\log( \, M \; / \; M^{}_0 \, )$", fontsize=18)
     ax.set_xlabel(r"$t \; / \; \tau^{}_{\rm KH}$", fontsize=18)
     ax.minorticks_on()
-    ax.tick_params(which='major', direction="in")
-    ax.tick_params(which='minor', direction="in")
+    ax.tick_params(which="major", direction="in")
+    ax.tick_params(which="minor", direction="in")
     ax.tick_params(labelsize=14)
 
     plt.savefig("mode_growth.png", dpi=300, bbox_inches="tight")
diff --git a/examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotSnapshots.py b/examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotSnapshots.py
index 391ab359585e75e7804c23d8d9c248ece25756a3..13b1f74a7b693a4b8590819e2f93c03a9f96869a 100644
--- a/examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotSnapshots.py
+++ b/examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotSnapshots.py
@@ -27,6 +27,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def make_axes():
     # Use Gridspec to set up figure
     n_gs_ax = 40
@@ -48,11 +49,33 @@ def make_axes():
     )
 
     ax_0 = plt.subplot(gs[:n_gs_ax, :n_gs_ax])
-    ax_1 = plt.subplot(gs[:n_gs_ax, n_gs_ax+n_gs_ax_gap:2*n_gs_ax+n_gs_ax_gap])
-    ax_2 = plt.subplot(gs[:n_gs_ax, 2*n_gs_ax+2*n_gs_ax_gap:3*n_gs_ax+2*n_gs_ax_gap])
+    ax_1 = plt.subplot(gs[:n_gs_ax, n_gs_ax + n_gs_ax_gap : 2 * n_gs_ax + n_gs_ax_gap])
+    ax_2 = plt.subplot(
+        gs[:n_gs_ax, 2 * n_gs_ax + 2 * n_gs_ax_gap : 3 * n_gs_ax + 2 * n_gs_ax_gap]
+    )
 
-    cax_0 = plt.subplot(gs[:int(0.5*n_gs_ax)-1, 3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap:3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap+n_gs_cbar])
-    cax_1 = plt.subplot(gs[int(0.5*n_gs_ax)+1:, 3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap:3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap+n_gs_cbar])
+    cax_0 = plt.subplot(
+        gs[
+            : int(0.5 * n_gs_ax) - 1,
+            3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
+    cax_1 = plt.subplot(
+        gs[
+            int(0.5 * n_gs_ax) + 1 :,
+            3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
 
     axs = [ax_0, ax_1, ax_2]
     caxs = [cax_0, cax_1]
@@ -67,8 +90,8 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
 
     with h5py.File(snap_file, "r") as f:
         # Units from file metadata to SI
-        m=float(f["Units"].attrs["Unit mass in cgs (U_M)"][0]) * 1e-3
-        l=float(f["Units"].attrs["Unit length in cgs (U_L)"][0]) * 1e-2
+        m = float(f["Units"].attrs["Unit mass in cgs (U_M)"][0]) * 1e-3
+        l = float(f["Units"].attrs["Unit length in cgs (U_L)"][0]) * 1e-2
 
         boxsize_l = f["Header"].attrs["BoxSize"][0] * l
         A1_x = f["/PartType0/Coordinates"][:, 0] * l
@@ -137,21 +160,20 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
     ax.set_ylim((0, boxsize_l))
 
 
-
 if __name__ == "__main__":
 
     # Set colormap
-    cmap1 = plt.get_cmap('YlOrRd')
+    cmap1 = plt.get_cmap("YlOrRd")
     mat_id1 = 402
     rho_min1 = 7950
     rho_max1 = 10050
-    norm1 = mpl.colors.Normalize(vmin=rho_min1,vmax=rho_max1)
+    norm1 = mpl.colors.Normalize(vmin=rho_min1, vmax=rho_max1)
 
-    cmap2 = plt.get_cmap('Blues_r')
+    cmap2 = plt.get_cmap("Blues_r")
     mat_id2 = 400
     rho_min2 = 4950
     rho_max2 = 5550
-    norm2 = mpl.colors.Normalize(vmin=rho_min2,vmax=rho_max2)
+    norm2 = mpl.colors.Normalize(vmin=rho_min2, vmax=rho_max2)
 
     # Generate axes
     axs, caxs = make_axes()
@@ -166,7 +188,14 @@ if __name__ == "__main__":
         time = times[i]
 
         plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2)
-        ax.text(0.5,-0.1,r"$t =\;$"+time+r"$\, \tau^{}_{\rm KH}$", horizontalalignment='center',size=18, transform=ax.transAxes)
+        ax.text(
+            0.5,
+            -0.1,
+            r"$t =\;$" + time + r"$\, \tau^{}_{\rm KH}$",
+            horizontalalignment="center",
+            size=18,
+            transform=ax.transAxes,
+        )
 
     # Colour bar
     sm1 = plt.cm.ScalarMappable(cmap=cmap1, norm=norm1)
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_3D/makeIC.py b/examples/Planetary/KelvinHelmholtz_IdealGas_3D/makeIC.py
index 3e6966579b7d6e0ff02acf90155e1db1469735b0..6169419a17e3be504a138bdeef39252f6a66afcb 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_3D/makeIC.py
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_3D/makeIC.py
@@ -1,22 +1,22 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
- #               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) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+#               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/>.
+#
+##############################################################################
 
 import h5py
 import numpy as np
@@ -24,22 +24,22 @@ import numpy as np
 # Generates a swift IC file for the Kelvin-Helmholtz test in a periodic box
 
 # Parameters
-N2_l     = 128    # Particles along one edge in the low-density region
-N2_depth = 18     # Particles in z direction in low-density region
-gamma = 5.0 / 3.0 # Gas adiabatic index
-matID1 = 0        # Central region material ID: ideal gas
-matID2 = 0        # Outskirts material ID: ideal gas
-P1    = 2.5       # Central region pressure
-P2    = 2.5       # Outskirts pressure
-rho1_approx =  2  # Central region density. Readjusted later
-rho2 =  1         # Outskirts density
-v1   = 0.5        # Central region velocity
-v2   = -0.5       # Outskirts velocity
-boxsize_l = 1     # size of simulation box in x and y dimension
-boxsize_depth = boxsize_l * N2_depth / N2_l # size of simulation box in z dimension
+N2_l = 128  # Particles along one edge in the low-density region
+N2_depth = 18  # Particles in z direction in low-density region
+gamma = 5.0 / 3.0  # Gas adiabatic index
+matID1 = 0  # Central region material ID: ideal gas
+matID2 = 0  # Outskirts material ID: ideal gas
+P1 = 2.5  # Central region pressure
+P2 = 2.5  # Outskirts pressure
+rho1_approx = 2  # Central region density. Readjusted later
+rho2 = 1  # Outskirts density
+v1 = 0.5  # Central region velocity
+v2 = -0.5  # Outskirts velocity
+boxsize_l = 1  # size of simulation box in x and y dimension
+boxsize_depth = boxsize_l * N2_depth / N2_l  # size of simulation box in z dimension
 mass = rho2 * (boxsize_l * boxsize_l * boxsize_depth) / (N2_l * N2_l * N2_depth)
 fileOutputName = "kelvin_helmholtz.hdf5"
-#---------------------------------------------------
+# ---------------------------------------------------
 
 # Start by calculating N1_l and rho1
 numPart2 = N2_l * N2_l * N2_depth
@@ -82,23 +82,27 @@ for i in range(N1_depth):
     for j in range(N1_l):
         for k in range(N1_l):
             index = i * N1_l**2 + j * N1_l + k
-            A2_coords1[index, 0] = (j / float(N1_l) + 1. / (2. * N1_l)) * boxsize_l
-            A2_coords1[index, 1] = (k / float(N1_l) + 1. / (2. * N1_l)) * boxsize_l
-            A2_coords1[index, 2] = (i / float(N1_depth) + 1. / (2. * N1_depth)) * boxsize_depth
+            A2_coords1[index, 0] = (j / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_l
+            A2_coords1[index, 1] = (k / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_l
+            A2_coords1[index, 2] = (
+                i / float(N1_depth) + 1.0 / (2.0 * N1_depth)
+            ) * boxsize_depth
 
 # Particles in the outskirts
 for i in range(N2_depth):
     for j in range(N2_l):
         for k in range(N2_l):
             index = i * N2_l**2 + j * N2_l + k
-            A2_coords2[index, 0] = (j / float(N2_l) + 1. / (2. * N2_l)) * boxsize_l
-            A2_coords2[index, 1] = (k / float(N2_l) + 1. / (2. * N2_l)) * boxsize_l
-            A2_coords2[index, 2] = (i / float(N2_depth) + 1. / (2. * N2_depth)) * boxsize_depth
+            A2_coords2[index, 0] = (j / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_l
+            A2_coords2[index, 1] = (k / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_l
+            A2_coords2[index, 2] = (
+                i / float(N2_depth) + 1.0 / (2.0 * N2_depth)
+            ) * boxsize_depth
 
 
 # Masks for the particles to be selected for the outer and inner regions
-mask1 = abs(A2_coords1[:,1] - 0.5 * boxsize_l) < 0.25 * boxsize_l
-mask2 = abs(A2_coords2[:,1] - 0.5 * boxsize_l) > 0.25 * boxsize_l
+mask1 = abs(A2_coords1[:, 1] - 0.5 * boxsize_l) < 0.25 * boxsize_l
+mask2 = abs(A2_coords2[:, 1] - 0.5 * boxsize_l) > 0.25 * boxsize_l
 
 # The positions of the particles are now selected
 # and the placement of the lattices are adjusted to give appropriate interfaces
@@ -111,16 +115,22 @@ pcl_separation_2 = np.cbrt(mass / rho2)
 boundary_separation = 0.5 * (pcl_separation_1 + pcl_separation_2)
 
 # Shift all the "inside" particles to get boundary_separation across the bottom interface
-min_y_inside = np.min(A2_coords_inside[:,1])
-max_y_outside_bot = np.max(A2_coords_outside[A2_coords_outside[:,1]-0.5 * boxsize_l < -0.25 * boxsize_l, 1])
+min_y_inside = np.min(A2_coords_inside[:, 1])
+max_y_outside_bot = np.max(
+    A2_coords_outside[A2_coords_outside[:, 1] - 0.5 * boxsize_l < -0.25 * boxsize_l, 1]
+)
 shift_distance_bot = boundary_separation - (min_y_inside - max_y_outside_bot)
-A2_coords_inside[:,1] += shift_distance_bot
+A2_coords_inside[:, 1] += shift_distance_bot
 
 # Shift the top section of the "outside" particles to get boundary_separation across the top interface
-max_y_inside = np.max(A2_coords_inside[:,1])
-min_y_outside_top = np.min(A2_coords_outside[A2_coords_outside[:,1]-0.5 * boxsize_l > 0.25 * boxsize_l, 1])
+max_y_inside = np.max(A2_coords_inside[:, 1])
+min_y_outside_top = np.min(
+    A2_coords_outside[A2_coords_outside[:, 1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1]
+)
 shift_distance_top = boundary_separation - (min_y_outside_top - max_y_inside)
-A2_coords_outside[A2_coords_outside[:,1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1] += shift_distance_top
+A2_coords_outside[
+    A2_coords_outside[:, 1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1
+] += shift_distance_top
 
 # Adjust box size in y direction based on the shifting of the lattices.
 new_box_y = boxsize_l + shift_distance_top
@@ -138,14 +148,16 @@ A1_ids = np.linspace(1, numPart, numPart)
 
 # Finally add the velocity perturbation
 vel_perturb_factor = 0.01 * (v1 - v2)
-A2_vel[:,1] = vel_perturb_factor * np.sin(2 * np.pi * A2_coords[:,0] / (0.5 * boxsize_l))
+A2_vel[:, 1] = vel_perturb_factor * np.sin(
+    2 * np.pi * A2_coords[:, 0] / (0.5 * boxsize_l)
+)
 
 # Write ICs to file
 with h5py.File(fileOutputName, "w") as f:
     # Header
     grp = f.create_group("/Header")
     grp.attrs["BoxSize"] = [boxsize_l, new_box_y, boxsize_depth]
-    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
@@ -154,29 +166,29 @@ with h5py.File(fileOutputName, "w") as f:
     grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
     grp.attrs["Dimension"] = 3
 
-    #Units
+    # Units
     grp = f.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 = f.create_group("/PartType0")
-    ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+    ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
     ds[()] = A2_coords
-    ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+    ds = grp.create_dataset("Velocities", (numPart, 3), "f")
     ds[()] = A2_vel
-    ds = grp.create_dataset('Masses', (numPart, 1), 'f')
-    ds[()] = A1_m.reshape((numPart,1))
-    ds = grp.create_dataset('Density', (numPart,1), 'f')
-    ds[()] = A1_rho.reshape((numPart,1))
-    ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-    ds[()] = A1_h.reshape((numPart,1))
-    ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-    ds[()] = A1_u.reshape((numPart,1))
-    ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
-    ds[()] = A1_ids.reshape((numPart,1))
-    ds = grp.create_dataset('MaterialIDs', (numPart,1), 'i')
-    ds[()] = A1_mat.reshape((numPart,1))
+    ds = grp.create_dataset("Masses", (numPart, 1), "f")
+    ds[()] = A1_m.reshape((numPart, 1))
+    ds = grp.create_dataset("Density", (numPart, 1), "f")
+    ds[()] = A1_rho.reshape((numPart, 1))
+    ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+    ds[()] = A1_h.reshape((numPart, 1))
+    ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+    ds[()] = A1_u.reshape((numPart, 1))
+    ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+    ds[()] = A1_ids.reshape((numPart, 1))
+    ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
+    ds[()] = A1_mat.reshape((numPart, 1))
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_3D/plotModeGrowth.py b/examples/Planetary/KelvinHelmholtz_IdealGas_3D/plotModeGrowth.py
index dab30dd07954f9b59e536be9a7b94511169bd2cb..1cfdf20fe5f85a20ee009466e3a1f70a0894eb63 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_3D/plotModeGrowth.py
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_3D/plotModeGrowth.py
@@ -28,6 +28,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def calculate_mode_growth(snap, vy_init_amp, wavelength):
 
     # Load data
@@ -40,46 +41,72 @@ def calculate_mode_growth(snap, vy_init_amp, wavelength):
         A1_h = f["/PartType0/SmoothingLengths"][:] / boxsize_l
 
     # Adjust positions
-    A2_pos[:,0] += 0.5
-    A2_pos[:,1] += 0.5
+    A2_pos[:, 0] += 0.5
+    A2_pos[:, 1] += 0.5
 
     # Masks to select the upper and lower halfs of the simulations
-    mask_up = A2_pos[:,1] >= 0.5
-    mask_down = A2_pos[:,1] < 0.5
+    mask_up = A2_pos[:, 1] >= 0.5
+    mask_down = A2_pos[:, 1] < 0.5
 
     # McNally et al. 2012 Eqn. 10
     s = np.empty(len(A1_h))
-    s[mask_down] = A2_vel[mask_down, 1] * A1_h[mask_down]**3 * np.sin(2 * np.pi * A2_pos[mask_down, 0] / wavelength) * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    s[mask_up] = A2_vel[mask_up, 1] * A1_h[mask_up]**3 * np.sin(2 * np.pi * A2_pos[mask_up, 0] / wavelength) * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    s[mask_down] = (
+        A2_vel[mask_down, 1]
+        * A1_h[mask_down] ** 3
+        * np.sin(2 * np.pi * A2_pos[mask_down, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
+    )
+    s[mask_up] = (
+        A2_vel[mask_up, 1]
+        * A1_h[mask_up] ** 3
+        * np.sin(2 * np.pi * A2_pos[mask_up, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    )
 
     # McNally et al. 2012 Eqn. 11
     c = np.empty(len(A1_h))
-    c[mask_down] = A2_vel[mask_down, 1] * A1_h[mask_down]**3 * np.cos(2 * np.pi * A2_pos[mask_down, 0] / wavelength) * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    c[mask_up] = A2_vel[mask_up, 1] * A1_h[mask_up]**3 * np.cos(2 * np.pi * A2_pos[mask_up, 0] / wavelength) * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    c[mask_down] = (
+        A2_vel[mask_down, 1]
+        * A1_h[mask_down] ** 3
+        * np.cos(2 * np.pi * A2_pos[mask_down, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
+    )
+    c[mask_up] = (
+        A2_vel[mask_up, 1]
+        * A1_h[mask_up] ** 3
+        * np.cos(2 * np.pi * A2_pos[mask_up, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    )
 
     # McNally et al. 2012 Eqn. 12
     d = np.empty(len(A1_h))
-    d[mask_down] = A1_h[mask_down]**3 * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    d[mask_up] = A1_h[mask_up]**3 * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    d[mask_down] = A1_h[mask_down] ** 3 * np.exp(
+        -2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength
+    )
+    d[mask_up] = A1_h[mask_up] ** 3 * np.exp(
+        -2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength
+    )
 
     # McNally et al. 2012 Eqn. 13
-    M = 2 * np.sqrt((np.sum(s)/ np.sum(d))**2 + (np.sum(c)/ np.sum(d))**2)
+    M = 2 * np.sqrt((np.sum(s) / np.sum(d)) ** 2 + (np.sum(c) / np.sum(d)) ** 2)
 
     return M
 
 
-
 if __name__ == "__main__":
 
     # Simulation paramerters for nomralisation of mode growth
-    vy_init_amp = 0.01 # Initial amplitude of y velocity perturbation in units of the boxsize per second
-    wavelength = 0.5   # wavelength of initial perturbation in units of the boxsize
+    vy_init_amp = 0.01  # Initial amplitude of y velocity perturbation in units of the boxsize per second
+    wavelength = 0.5  # wavelength of initial perturbation in units of the boxsize
 
     # Use Gridspec to set up figure
     n_gs_ax = 40
     ax_len = 5
     fig = plt.figure(figsize=(9, 6))
-    gs = mpl.gridspec.GridSpec(nrows=40, ncols=60,)
+    gs = mpl.gridspec.GridSpec(
+        nrows=40,
+        ncols=60,
+    )
     ax = plt.subplot(gs[:, :])
 
     # Snapshots and corresponding times
@@ -98,8 +125,8 @@ if __name__ == "__main__":
     ax.set_ylabel(r"$\log( \, M \; / \; M^{}_0 \, )$", fontsize=18)
     ax.set_xlabel(r"$t \; / \; \tau^{}_{\rm KH}$", fontsize=18)
     ax.minorticks_on()
-    ax.tick_params(which='major', direction="in")
-    ax.tick_params(which='minor', direction="in")
+    ax.tick_params(which="major", direction="in")
+    ax.tick_params(which="minor", direction="in")
     ax.tick_params(labelsize=14)
 
     plt.savefig("mode_growth.png", dpi=300, bbox_inches="tight")
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_3D/plotSnapshots.py b/examples/Planetary/KelvinHelmholtz_IdealGas_3D/plotSnapshots.py
index 85df472ec90b7f51f48d0c67cb22abdb5d07bf7d..0fbcade7d99150c4cc6b298c4de40df2e2e203fd 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_3D/plotSnapshots.py
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_3D/plotSnapshots.py
@@ -27,6 +27,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def make_axes():
     # Use Gridspec to set up figure
     n_gs_ax = 40
@@ -48,9 +49,21 @@ def make_axes():
     )
 
     ax_0 = plt.subplot(gs[:n_gs_ax, :n_gs_ax])
-    ax_1 = plt.subplot(gs[:n_gs_ax, n_gs_ax+n_gs_ax_gap:2*n_gs_ax+n_gs_ax_gap])
-    ax_2 = plt.subplot(gs[:n_gs_ax, 2*n_gs_ax+2*n_gs_ax_gap:3*n_gs_ax+2*n_gs_ax_gap])
-    cax = plt.subplot(gs[:n_gs_ax, 3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap:3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap+n_gs_cbar])
+    ax_1 = plt.subplot(gs[:n_gs_ax, n_gs_ax + n_gs_ax_gap : 2 * n_gs_ax + n_gs_ax_gap])
+    ax_2 = plt.subplot(
+        gs[:n_gs_ax, 2 * n_gs_ax + 2 * n_gs_ax_gap : 3 * n_gs_ax + 2 * n_gs_ax_gap]
+    )
+    cax = plt.subplot(
+        gs[
+            :n_gs_ax,
+            3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
 
     axs = [ax_0, ax_1, ax_2]
 
@@ -113,7 +126,6 @@ def plot_kh(ax, snap, cmap, norm):
     ax.set_ylim((0, boxsize_l))
 
 
-
 if __name__ == "__main__":
 
     # Set colormap
@@ -134,7 +146,14 @@ if __name__ == "__main__":
         time = times[i]
 
         plot_kh(ax, snap, cmap, norm)
-        ax.text(0.5,-0.1,r"$t =\;$"+time+r"$\, \tau^{}_{\rm KH}$", horizontalalignment='center',size=18, transform=ax.transAxes)
+        ax.text(
+            0.5,
+            -0.1,
+            r"$t =\;$" + time + r"$\, \tau^{}_{\rm KH}$",
+            horizontalalignment="center",
+            size=18,
+            transform=ax.transAxes,
+        )
 
     # Colour bar
     sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/makeIC.py b/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/makeIC.py
index 2a28d194e7fa315d3f181e9edf6fa1f3607b50af..9f14d9f76cf50184c31585396d30322fb8d60037 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/makeIC.py
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/makeIC.py
@@ -1,22 +1,22 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
- #               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) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+#               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/>.
+#
+##############################################################################
 
 import h5py
 import numpy as np
@@ -24,22 +24,22 @@ import numpy as np
 # Generates a swift IC file for the Kelvin-Helmholtz test in a periodic box
 
 # Parameters
-N2_l     = 256    # Particles along one edge in the low-density region
-N2_depth = 18     # Particles in z direction in low-density region
-gamma = 5.0 / 3.0 # Gas adiabatic index
-matID1 = 0        # Central region material ID: ideal gas
-matID2 = 0        # Outskirts material ID: ideal gas
-P1    = 2.5       # Central region pressure
-P2    = 2.5       # Outskirts pressure
-rho1_approx =  11 # Central region density. Readjusted later
-rho2 =  1         # Outskirts density
-v1   = 0.5        # Central region velocity
-v2   = -0.5       # Outskirts velocity
-boxsize_l = 1     # size of simulation box in x and y dimension
-boxsize_depth = boxsize_l * N2_depth / N2_l # size of simulation box in z dimension
+N2_l = 256  # Particles along one edge in the low-density region
+N2_depth = 18  # Particles in z direction in low-density region
+gamma = 5.0 / 3.0  # Gas adiabatic index
+matID1 = 0  # Central region material ID: ideal gas
+matID2 = 0  # Outskirts material ID: ideal gas
+P1 = 2.5  # Central region pressure
+P2 = 2.5  # Outskirts pressure
+rho1_approx = 11  # Central region density. Readjusted later
+rho2 = 1  # Outskirts density
+v1 = 0.5  # Central region velocity
+v2 = -0.5  # Outskirts velocity
+boxsize_l = 1  # size of simulation box in x and y dimension
+boxsize_depth = boxsize_l * N2_depth / N2_l  # size of simulation box in z dimension
 mass = rho2 * (boxsize_l * boxsize_l * boxsize_depth) / (N2_l * N2_l * N2_depth)
 fileOutputName = "kelvin_helmholtz.hdf5"
-#---------------------------------------------------
+# ---------------------------------------------------
 
 # Start by calculating N1_l and rho1
 numPart2 = N2_l * N2_l * N2_depth
@@ -82,23 +82,27 @@ for i in range(N1_depth):
     for j in range(N1_l):
         for k in range(N1_l):
             index = i * N1_l**2 + j * N1_l + k
-            A2_coords1[index, 0] = (j / float(N1_l) + 1. / (2. * N1_l)) * boxsize_l
-            A2_coords1[index, 1] = (k / float(N1_l) + 1. / (2. * N1_l)) * boxsize_l
-            A2_coords1[index, 2] = (i / float(N1_depth) + 1. / (2. * N1_depth)) * boxsize_depth
+            A2_coords1[index, 0] = (j / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_l
+            A2_coords1[index, 1] = (k / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_l
+            A2_coords1[index, 2] = (
+                i / float(N1_depth) + 1.0 / (2.0 * N1_depth)
+            ) * boxsize_depth
 
 # Particles in the outskirts
 for i in range(N2_depth):
     for j in range(N2_l):
         for k in range(N2_l):
             index = i * N2_l**2 + j * N2_l + k
-            A2_coords2[index, 0] = (j / float(N2_l) + 1. / (2. * N2_l)) * boxsize_l
-            A2_coords2[index, 1] = (k / float(N2_l) + 1. / (2. * N2_l)) * boxsize_l
-            A2_coords2[index, 2] = (i / float(N2_depth) + 1. / (2. * N2_depth)) * boxsize_depth
+            A2_coords2[index, 0] = (j / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_l
+            A2_coords2[index, 1] = (k / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_l
+            A2_coords2[index, 2] = (
+                i / float(N2_depth) + 1.0 / (2.0 * N2_depth)
+            ) * boxsize_depth
 
 
 # Masks for the particles to be selected for the outer and inner regions
-mask1 = abs(A2_coords1[:,1] - 0.5 * boxsize_l) < 0.25 * boxsize_l
-mask2 = abs(A2_coords2[:,1] - 0.5 * boxsize_l) > 0.25 * boxsize_l
+mask1 = abs(A2_coords1[:, 1] - 0.5 * boxsize_l) < 0.25 * boxsize_l
+mask2 = abs(A2_coords2[:, 1] - 0.5 * boxsize_l) > 0.25 * boxsize_l
 
 # The positions of the particles are now selected
 # and the placement of the lattices are adjusted to give appropriate interfaces
@@ -111,16 +115,22 @@ pcl_separation_2 = np.cbrt(mass / rho2)
 boundary_separation = 0.5 * (pcl_separation_1 + pcl_separation_2)
 
 # Shift all the "inside" particles to get boundary_separation across the bottom interface
-min_y_inside = np.min(A2_coords_inside[:,1])
-max_y_outside_bot = np.max(A2_coords_outside[A2_coords_outside[:,1]-0.5 * boxsize_l < -0.25 * boxsize_l, 1])
+min_y_inside = np.min(A2_coords_inside[:, 1])
+max_y_outside_bot = np.max(
+    A2_coords_outside[A2_coords_outside[:, 1] - 0.5 * boxsize_l < -0.25 * boxsize_l, 1]
+)
 shift_distance_bot = boundary_separation - (min_y_inside - max_y_outside_bot)
-A2_coords_inside[:,1] += shift_distance_bot
+A2_coords_inside[:, 1] += shift_distance_bot
 
 # Shift the top section of the "outside" particles to get boundary_separation across the top interface
-max_y_inside = np.max(A2_coords_inside[:,1])
-min_y_outside_top = np.min(A2_coords_outside[A2_coords_outside[:,1]-0.5 * boxsize_l > 0.25 * boxsize_l, 1])
+max_y_inside = np.max(A2_coords_inside[:, 1])
+min_y_outside_top = np.min(
+    A2_coords_outside[A2_coords_outside[:, 1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1]
+)
 shift_distance_top = boundary_separation - (min_y_outside_top - max_y_inside)
-A2_coords_outside[A2_coords_outside[:,1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1] += shift_distance_top
+A2_coords_outside[
+    A2_coords_outside[:, 1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1
+] += shift_distance_top
 
 # Adjust box size in y direction based on the shifting of the lattices.
 new_box_y = boxsize_l + shift_distance_top
@@ -138,14 +148,16 @@ A1_ids = np.linspace(1, numPart, numPart)
 
 # Finally add the velocity perturbation
 vel_perturb_factor = 0.01 * (v1 - v2)
-A2_vel[:,1] = vel_perturb_factor * np.sin(2 * np.pi * A2_coords[:,0] / (0.125 * boxsize_l))
+A2_vel[:, 1] = vel_perturb_factor * np.sin(
+    2 * np.pi * A2_coords[:, 0] / (0.125 * boxsize_l)
+)
 
 # Write ICs to file
 with h5py.File(fileOutputName, "w") as f:
     # Header
     grp = f.create_group("/Header")
     grp.attrs["BoxSize"] = [boxsize_l, new_box_y, boxsize_depth]
-    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
@@ -154,29 +166,29 @@ with h5py.File(fileOutputName, "w") as f:
     grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
     grp.attrs["Dimension"] = 3
 
-    #Units
+    # Units
     grp = f.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 = f.create_group("/PartType0")
-    ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+    ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
     ds[()] = A2_coords
-    ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+    ds = grp.create_dataset("Velocities", (numPart, 3), "f")
     ds[()] = A2_vel
-    ds = grp.create_dataset('Masses', (numPart, 1), 'f')
-    ds[()] = A1_m.reshape((numPart,1))
-    ds = grp.create_dataset('Density', (numPart,1), 'f')
-    ds[()] = A1_rho.reshape((numPart,1))
-    ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-    ds[()] = A1_h.reshape((numPart,1))
-    ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-    ds[()] = A1_u.reshape((numPart,1))
-    ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
-    ds[()] = A1_ids.reshape((numPart,1))
-    ds = grp.create_dataset('MaterialIDs', (numPart,1), 'i')
-    ds[()] = A1_mat.reshape((numPart,1))
+    ds = grp.create_dataset("Masses", (numPart, 1), "f")
+    ds[()] = A1_m.reshape((numPart, 1))
+    ds = grp.create_dataset("Density", (numPart, 1), "f")
+    ds[()] = A1_rho.reshape((numPart, 1))
+    ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+    ds[()] = A1_h.reshape((numPart, 1))
+    ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+    ds[()] = A1_u.reshape((numPart, 1))
+    ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+    ds[()] = A1_ids.reshape((numPart, 1))
+    ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
+    ds[()] = A1_mat.reshape((numPart, 1))
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/plotModeGrowth.py b/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/plotModeGrowth.py
index 10c8b8ca187cf590075bc14920bb3fc900ef2062..14ec1920af563ba5fb70a729efc059d9e7deba9b 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/plotModeGrowth.py
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/plotModeGrowth.py
@@ -28,6 +28,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def calculate_mode_growth(snap, vy_init_amp, wavelength):
 
     # Load data
@@ -40,46 +41,72 @@ def calculate_mode_growth(snap, vy_init_amp, wavelength):
         A1_h = f["/PartType0/SmoothingLengths"][:] / boxsize_l
 
     # Adjust positions
-    A2_pos[:,0] += 0.5
-    A2_pos[:,1] += 0.5
+    A2_pos[:, 0] += 0.5
+    A2_pos[:, 1] += 0.5
 
     # Masks to select the upper and lower halfs of the simulations
-    mask_up = A2_pos[:,1] >= 0.5
-    mask_down = A2_pos[:,1] < 0.5
+    mask_up = A2_pos[:, 1] >= 0.5
+    mask_down = A2_pos[:, 1] < 0.5
 
     # McNally et al. 2012 Eqn. 10
     s = np.empty(len(A1_h))
-    s[mask_down] = A2_vel[mask_down, 1] * A1_h[mask_down]**3 * np.sin(2 * np.pi * A2_pos[mask_down, 0] / wavelength) * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    s[mask_up] = A2_vel[mask_up, 1] * A1_h[mask_up]**3 * np.sin(2 * np.pi * A2_pos[mask_up, 0] / wavelength) * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    s[mask_down] = (
+        A2_vel[mask_down, 1]
+        * A1_h[mask_down] ** 3
+        * np.sin(2 * np.pi * A2_pos[mask_down, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
+    )
+    s[mask_up] = (
+        A2_vel[mask_up, 1]
+        * A1_h[mask_up] ** 3
+        * np.sin(2 * np.pi * A2_pos[mask_up, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    )
 
     # McNally et al. 2012 Eqn. 11
     c = np.empty(len(A1_h))
-    c[mask_down] = A2_vel[mask_down, 1] * A1_h[mask_down]**3 * np.cos(2 * np.pi * A2_pos[mask_down, 0] / wavelength) * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    c[mask_up] = A2_vel[mask_up, 1] * A1_h[mask_up]**3 * np.cos(2 * np.pi * A2_pos[mask_up, 0] / wavelength) * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    c[mask_down] = (
+        A2_vel[mask_down, 1]
+        * A1_h[mask_down] ** 3
+        * np.cos(2 * np.pi * A2_pos[mask_down, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
+    )
+    c[mask_up] = (
+        A2_vel[mask_up, 1]
+        * A1_h[mask_up] ** 3
+        * np.cos(2 * np.pi * A2_pos[mask_up, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    )
 
     # McNally et al. 2012 Eqn. 12
     d = np.empty(len(A1_h))
-    d[mask_down] = A1_h[mask_down]**3 * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    d[mask_up] = A1_h[mask_up]**3 * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    d[mask_down] = A1_h[mask_down] ** 3 * np.exp(
+        -2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength
+    )
+    d[mask_up] = A1_h[mask_up] ** 3 * np.exp(
+        -2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength
+    )
 
     # McNally et al. 2012 Eqn. 13
-    M = 2 * np.sqrt((np.sum(s)/ np.sum(d))**2 + (np.sum(c)/ np.sum(d))**2)
+    M = 2 * np.sqrt((np.sum(s) / np.sum(d)) ** 2 + (np.sum(c) / np.sum(d)) ** 2)
 
     return M
 
 
-
 if __name__ == "__main__":
 
     # Simulation paramerters for nomralisation of mode growth
-    vy_init_amp = 0.01 # Initial amplitude of y velocity perturbation in units of the boxsize per second
-    wavelength = 0.125   # wavelength of initial perturbation in units of the boxsize
+    vy_init_amp = 0.01  # Initial amplitude of y velocity perturbation in units of the boxsize per second
+    wavelength = 0.125  # wavelength of initial perturbation in units of the boxsize
 
     # Use Gridspec to set up figure
     n_gs_ax = 40
     ax_len = 5
     fig = plt.figure(figsize=(9, 6))
-    gs = mpl.gridspec.GridSpec(nrows=40, ncols=60,)
+    gs = mpl.gridspec.GridSpec(
+        nrows=40,
+        ncols=60,
+    )
     ax = plt.subplot(gs[:, :])
 
     # Snapshots and corresponding times
@@ -98,8 +125,8 @@ if __name__ == "__main__":
     ax.set_ylabel(r"$\log( \, M \; / \; M^{}_0 \, )$", fontsize=18)
     ax.set_xlabel(r"$t \; / \; \tau^{}_{\rm KH}$", fontsize=18)
     ax.minorticks_on()
-    ax.tick_params(which='major', direction="in")
-    ax.tick_params(which='minor', direction="in")
+    ax.tick_params(which="major", direction="in")
+    ax.tick_params(which="minor", direction="in")
     ax.tick_params(labelsize=14)
 
     plt.savefig("mode_growth.png", dpi=300, bbox_inches="tight")
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/plotSnapshots.py b/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/plotSnapshots.py
index 2f5d1404d765b9f3475c327aae5f643b0ab097a4..52778f73fab84f8627d675653cc878b30e4440df 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/plotSnapshots.py
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_LargeDensityContrast_3D/plotSnapshots.py
@@ -27,6 +27,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def make_axes():
     # Use Gridspec to set up figure
     n_gs_ax = 40
@@ -48,9 +49,21 @@ def make_axes():
     )
 
     ax_0 = plt.subplot(gs[:n_gs_ax, :n_gs_ax])
-    ax_1 = plt.subplot(gs[:n_gs_ax, n_gs_ax+n_gs_ax_gap:2*n_gs_ax+n_gs_ax_gap])
-    ax_2 = plt.subplot(gs[:n_gs_ax, 2*n_gs_ax+2*n_gs_ax_gap:3*n_gs_ax+2*n_gs_ax_gap])
-    cax = plt.subplot(gs[:n_gs_ax, 3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap:3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap+n_gs_cbar])
+    ax_1 = plt.subplot(gs[:n_gs_ax, n_gs_ax + n_gs_ax_gap : 2 * n_gs_ax + n_gs_ax_gap])
+    ax_2 = plt.subplot(
+        gs[:n_gs_ax, 2 * n_gs_ax + 2 * n_gs_ax_gap : 3 * n_gs_ax + 2 * n_gs_ax_gap]
+    )
+    cax = plt.subplot(
+        gs[
+            :n_gs_ax,
+            3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
 
     axs = [ax_0, ax_1, ax_2]
 
@@ -113,13 +126,12 @@ def plot_kh(ax, snap, cmap, norm):
     ax.set_ylim((0, boxsize_l))
 
 
-
 if __name__ == "__main__":
 
     # Set colormap
     cmap = plt.get_cmap("Spectral_r")
     vmin, vmax = 0.95, 2.05
-    norm = mpl.colors.LogNorm(vmin=0.8,vmax=12)
+    norm = mpl.colors.LogNorm(vmin=0.8, vmax=12)
 
     # Generate axes
     axs, cax = make_axes()
@@ -134,7 +146,14 @@ if __name__ == "__main__":
         time = times[i]
 
         plot_kh(ax, snap, cmap, norm)
-        ax.text(0.5,-0.1,r"$t =\;$"+time+r"$\, \tau^{}_{\rm KH}$", horizontalalignment='center',size=18, transform=ax.transAxes)
+        ax.text(
+            0.5,
+            -0.1,
+            r"$t =\;$" + time + r"$\, \tau^{}_{\rm KH}$",
+            horizontalalignment="center",
+            size=18,
+            transform=ax.transAxes,
+        )
 
     # Colour bar
     sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/makeIC.py b/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/makeIC.py
index d01b0a257a52f847f732ea336ad1e221cca3a11a..ca0026d3c7d26eb55af39f98f690bfb62d7f9a6b 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/makeIC.py
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/makeIC.py
@@ -1,22 +1,22 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
- #               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) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+#               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/>.
+#
+##############################################################################
 
 import h5py
 import numpy as np
@@ -24,24 +24,24 @@ import numpy as np
 # Generates a swift IC file for the Kelvin-Helmholtz test in a periodic box
 
 # Parameters
-N_l     = 128    # Particles along one edge in the low-density region
-N_depth = 18     # Particles in z direction in low-density region
-gamma = 5.0 / 3.0 # Gas adiabatic index
-P1    = 2.5       # Central region pressure
-P2    = 2.5       # Outskirts pressure
-rho1 =  2         # Central region density
-rho2 =  1         # Outskirts density
-v1   = 0.5        # Central region velocity
-v2   = -0.5       # Outskirts velocity
-boxsize_l = 1     # size of simulation box in x and y dimension
-boxsize_depth = boxsize_l * N_depth / N_l # size of simulation box in z dimension
+N_l = 128  # Particles along one edge in the low-density region
+N_depth = 18  # Particles in z direction in low-density region
+gamma = 5.0 / 3.0  # Gas adiabatic index
+P1 = 2.5  # Central region pressure
+P2 = 2.5  # Outskirts pressure
+rho1 = 2  # Central region density
+rho2 = 1  # Outskirts density
+v1 = 0.5  # Central region velocity
+v2 = -0.5  # Outskirts velocity
+boxsize_l = 1  # size of simulation box in x and y dimension
+boxsize_depth = boxsize_l * N_depth / N_l  # size of simulation box in z dimension
 fileOutputName = "kelvin_helmholtz.hdf5"
 
 # Parameters for smoothing of interfaces
-vm    = (v2 - v1) / 2
+vm = (v2 - v1) / 2
 rhom = (rho2 - rho1) / 2
 delta = 0.025
-#---------------------------------------------------
+# ---------------------------------------------------
 
 numPart = N_l * N_l * N_depth
 
@@ -62,48 +62,50 @@ for i in range(N_depth):
         for k in range(N_l):
             index = i * N_l**2 + j * N_l + k
 
-            x = (j / float(N_l) + 1. / (2. * N_l)) * boxsize_l
-            y = (k / float(N_l) + 1. / (2. * N_l)) * boxsize_l
-            z = (i / float(N_depth) + 1. / (2. * N_depth)) * boxsize_depth
+            x = (j / float(N_l) + 1.0 / (2.0 * N_l)) * boxsize_l
+            y = (k / float(N_l) + 1.0 / (2.0 * N_l)) * boxsize_l
+            z = (i / float(N_depth) + 1.0 / (2.0 * N_depth)) * boxsize_depth
 
             A2_coords[index, 0] = x
             A2_coords[index, 1] = y
             A2_coords[index, 2] = z
 
-            if 0.0<=y<=0.25:
-                A1_rho[index] = rho2 - rhom * np.exp((y-0.25)/delta)
-                A2_vel[index, 0] = v2 - vm * np.exp((y-0.25)/delta)
+            if 0.0 <= y <= 0.25:
+                A1_rho[index] = rho2 - rhom * np.exp((y - 0.25) / delta)
+                A2_vel[index, 0] = v2 - vm * np.exp((y - 0.25) / delta)
                 A1_m[index] = A1_rho[index] / N_l**3
-                A1_u[index] = P2 / (A1_rho[index] * (gamma-1.))
+                A1_u[index] = P2 / (A1_rho[index] * (gamma - 1.0))
 
-            elif 0.25<=y<=0.5:
-                A1_rho[index] = rho1 + rhom * np.exp((0.25-y)/delta)
-                A2_vel[index, 0] = v1 + vm * np.exp((0.25-y)/delta)
+            elif 0.25 <= y <= 0.5:
+                A1_rho[index] = rho1 + rhom * np.exp((0.25 - y) / delta)
+                A2_vel[index, 0] = v1 + vm * np.exp((0.25 - y) / delta)
                 A1_m[index] = A1_rho[index] / N_l**3
-                A1_u[index] = P1 / (A1_rho[index] * (gamma-1.))
+                A1_u[index] = P1 / (A1_rho[index] * (gamma - 1.0))
 
-            elif 0.5<=y<=0.75:
-                A1_rho[index] = rho1 + rhom * np.exp((y-0.75)/delta)
-                A2_vel[index, 0] = v1 + vm * np.exp((y-0.75)/delta)
+            elif 0.5 <= y <= 0.75:
+                A1_rho[index] = rho1 + rhom * np.exp((y - 0.75) / delta)
+                A2_vel[index, 0] = v1 + vm * np.exp((y - 0.75) / delta)
                 A1_m[index] = A1_rho[index] / N_l**3
-                A1_u[index] = P1 / (A1_rho[index] * (gamma-1.))
+                A1_u[index] = P1 / (A1_rho[index] * (gamma - 1.0))
 
-            elif 0.75<=y<=1:
-                A1_rho[index] = rho2 - rhom * np.exp((0.75-y)/delta)
-                A2_vel[index, 0] = v2 - vm * np.exp((0.75-y)/delta)
+            elif 0.75 <= y <= 1:
+                A1_rho[index] = rho2 - rhom * np.exp((0.75 - y) / delta)
+                A2_vel[index, 0] = v2 - vm * np.exp((0.75 - y) / delta)
                 A1_m[index] = A1_rho[index] / N_l**3
-                A1_u[index] = P2 / (A1_rho[index] * (gamma-1.))
+                A1_u[index] = P2 / (A1_rho[index] * (gamma - 1.0))
 
 # Finally add the velocity perturbation
 vel_perturb_factor = 0.01 * (v1 - v2)
-A2_vel[:,1] = vel_perturb_factor * np.sin(2 * np.pi * A2_coords[:,0] / (0.5 * boxsize_l))
+A2_vel[:, 1] = vel_perturb_factor * np.sin(
+    2 * np.pi * A2_coords[:, 0] / (0.5 * boxsize_l)
+)
 
 # Write ICs to file
 with h5py.File(fileOutputName, "w") as f:
     # Header
     grp = f.create_group("/Header")
     grp.attrs["BoxSize"] = [boxsize_l, boxsize_l, boxsize_depth]
-    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
@@ -112,29 +114,29 @@ with h5py.File(fileOutputName, "w") as f:
     grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
     grp.attrs["Dimension"] = 3
 
-    #Units
+    # Units
     grp = f.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 = f.create_group("/PartType0")
-    ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+    ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
     ds[()] = A2_coords
-    ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+    ds = grp.create_dataset("Velocities", (numPart, 3), "f")
     ds[()] = A2_vel
-    ds = grp.create_dataset('Masses', (numPart, 1), 'f')
-    ds[()] = A1_m.reshape((numPart,1))
-    ds = grp.create_dataset('Density', (numPart,1), 'f')
-    ds[()] = A1_rho.reshape((numPart,1))
-    ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-    ds[()] = A1_h.reshape((numPart,1))
-    ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-    ds[()] = A1_u.reshape((numPart,1))
-    ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
-    ds[()] = A1_ids.reshape((numPart,1))
-    ds = grp.create_dataset('MaterialIDs', (numPart,1), 'i')
-    ds[()] = A1_mat.reshape((numPart,1))
+    ds = grp.create_dataset("Masses", (numPart, 1), "f")
+    ds[()] = A1_m.reshape((numPart, 1))
+    ds = grp.create_dataset("Density", (numPart, 1), "f")
+    ds[()] = A1_rho.reshape((numPart, 1))
+    ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+    ds[()] = A1_h.reshape((numPart, 1))
+    ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+    ds[()] = A1_u.reshape((numPart, 1))
+    ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+    ds[()] = A1_ids.reshape((numPart, 1))
+    ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
+    ds[()] = A1_mat.reshape((numPart, 1))
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/plotModeGrowth.py b/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/plotModeGrowth.py
index fde18517b8f2c945567d5cceb087880f815779df..59a4d6bcb85369979c35d924af55247f8693da41 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/plotModeGrowth.py
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/plotModeGrowth.py
@@ -28,6 +28,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def calculate_mode_growth(snap, vy_init_amp, wavelength):
 
     # Load data
@@ -40,46 +41,72 @@ def calculate_mode_growth(snap, vy_init_amp, wavelength):
         A1_h = f["/PartType0/SmoothingLengths"][:] / boxsize_l
 
     # Adjust positions
-    A2_pos[:,0] += 0.5
-    A2_pos[:,1] += 0.5
+    A2_pos[:, 0] += 0.5
+    A2_pos[:, 1] += 0.5
 
     # Masks to select the upper and lower halfs of the simulations
-    mask_up = A2_pos[:,1] >= 0.5
-    mask_down = A2_pos[:,1] < 0.5
+    mask_up = A2_pos[:, 1] >= 0.5
+    mask_down = A2_pos[:, 1] < 0.5
 
     # McNally et al. 2012 Eqn. 10
     s = np.empty(len(A1_h))
-    s[mask_down] = A2_vel[mask_down, 1] * A1_h[mask_down]**3 * np.sin(2 * np.pi * A2_pos[mask_down, 0] / wavelength) * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    s[mask_up] = A2_vel[mask_up, 1] * A1_h[mask_up]**3 * np.sin(2 * np.pi * A2_pos[mask_up, 0] / wavelength) * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    s[mask_down] = (
+        A2_vel[mask_down, 1]
+        * A1_h[mask_down] ** 3
+        * np.sin(2 * np.pi * A2_pos[mask_down, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
+    )
+    s[mask_up] = (
+        A2_vel[mask_up, 1]
+        * A1_h[mask_up] ** 3
+        * np.sin(2 * np.pi * A2_pos[mask_up, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    )
 
     # McNally et al. 2012 Eqn. 11
     c = np.empty(len(A1_h))
-    c[mask_down] = A2_vel[mask_down, 1] * A1_h[mask_down]**3 * np.cos(2 * np.pi * A2_pos[mask_down, 0] / wavelength) * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    c[mask_up] = A2_vel[mask_up, 1] * A1_h[mask_up]**3 * np.cos(2 * np.pi * A2_pos[mask_up, 0] / wavelength) * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    c[mask_down] = (
+        A2_vel[mask_down, 1]
+        * A1_h[mask_down] ** 3
+        * np.cos(2 * np.pi * A2_pos[mask_down, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
+    )
+    c[mask_up] = (
+        A2_vel[mask_up, 1]
+        * A1_h[mask_up] ** 3
+        * np.cos(2 * np.pi * A2_pos[mask_up, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    )
 
     # McNally et al. 2012 Eqn. 12
     d = np.empty(len(A1_h))
-    d[mask_down] = A1_h[mask_down]**3 * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    d[mask_up] = A1_h[mask_up]**3 * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    d[mask_down] = A1_h[mask_down] ** 3 * np.exp(
+        -2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength
+    )
+    d[mask_up] = A1_h[mask_up] ** 3 * np.exp(
+        -2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength
+    )
 
     # McNally et al. 2012 Eqn. 13
-    M = 2 * np.sqrt((np.sum(s)/ np.sum(d))**2 + (np.sum(c)/ np.sum(d))**2)
+    M = 2 * np.sqrt((np.sum(s) / np.sum(d)) ** 2 + (np.sum(c) / np.sum(d)) ** 2)
 
     return M
 
 
-
 if __name__ == "__main__":
 
     # Simulation paramerters for nomralisation of mode growth
-    vy_init_amp = 0.01 # Initial amplitude of y velocity perturbation in units of the boxsize per second
-    wavelength = 0.5   # wavelength of initial perturbation in units of the boxsize
+    vy_init_amp = 0.01  # Initial amplitude of y velocity perturbation in units of the boxsize per second
+    wavelength = 0.5  # wavelength of initial perturbation in units of the boxsize
 
     # Use Gridspec to set up figure
     n_gs_ax = 40
     ax_len = 5
     fig = plt.figure(figsize=(9, 6))
-    gs = mpl.gridspec.GridSpec(nrows=40, ncols=60,)
+    gs = mpl.gridspec.GridSpec(
+        nrows=40,
+        ncols=60,
+    )
     ax = plt.subplot(gs[:, :])
 
     # Snapshots and corresponding times
@@ -98,8 +125,8 @@ if __name__ == "__main__":
     ax.set_ylabel(r"$\log( \, M \; / \; M^{}_0 \, )$", fontsize=18)
     ax.set_xlabel(r"$t \; / \; \tau^{}_{\rm KH}$", fontsize=18)
     ax.minorticks_on()
-    ax.tick_params(which='major', direction="in")
-    ax.tick_params(which='minor', direction="in")
+    ax.tick_params(which="major", direction="in")
+    ax.tick_params(which="minor", direction="in")
     ax.tick_params(labelsize=14)
 
     plt.savefig("mode_growth.png", dpi=300, bbox_inches="tight")
diff --git a/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/plotSnapshots.py b/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/plotSnapshots.py
index 0c4394ecef7bd10626019076c176f636eba565e1..3ffbb3d6b38ed3d6a5b6eee5b49255d4fed20c58 100644
--- a/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/plotSnapshots.py
+++ b/examples/Planetary/KelvinHelmholtz_IdealGas_Smooth_3D/plotSnapshots.py
@@ -27,6 +27,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def make_axes():
     # Use Gridspec to set up figure
     n_gs_ax = 40
@@ -48,9 +49,21 @@ def make_axes():
     )
 
     ax_0 = plt.subplot(gs[:n_gs_ax, :n_gs_ax])
-    ax_1 = plt.subplot(gs[:n_gs_ax, n_gs_ax+n_gs_ax_gap:2*n_gs_ax+n_gs_ax_gap])
-    ax_2 = plt.subplot(gs[:n_gs_ax, 2*n_gs_ax+2*n_gs_ax_gap:3*n_gs_ax+2*n_gs_ax_gap])
-    cax = plt.subplot(gs[:n_gs_ax, 3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap:3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap+n_gs_cbar])
+    ax_1 = plt.subplot(gs[:n_gs_ax, n_gs_ax + n_gs_ax_gap : 2 * n_gs_ax + n_gs_ax_gap])
+    ax_2 = plt.subplot(
+        gs[:n_gs_ax, 2 * n_gs_ax + 2 * n_gs_ax_gap : 3 * n_gs_ax + 2 * n_gs_ax_gap]
+    )
+    cax = plt.subplot(
+        gs[
+            :n_gs_ax,
+            3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
 
     axs = [ax_0, ax_1, ax_2]
 
@@ -113,7 +126,6 @@ def plot_kh(ax, snap, cmap, norm):
     ax.set_ylim((0, boxsize_l))
 
 
-
 if __name__ == "__main__":
 
     # Set colormap
@@ -134,7 +146,14 @@ if __name__ == "__main__":
         time = times[i]
 
         plot_kh(ax, snap, cmap, norm)
-        ax.text(0.5,-0.1,r"$t =\;$"+time+r"$\, \tau^{}_{\rm KH}$", horizontalalignment='center',size=18, transform=ax.transAxes)
+        ax.text(
+            0.5,
+            -0.1,
+            r"$t =\;$" + time + r"$\, \tau^{}_{\rm KH}$",
+            horizontalalignment="center",
+            size=18,
+            transform=ax.transAxes,
+        )
 
     # Colour bar
     sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
diff --git a/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/makeIC.py b/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/makeIC.py
index 71e114f3b341170a78ed743ed1abef3635554676..721d04a3ade81516491d9cff3e2f7d27a4ab6fa7 100644
--- a/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/makeIC.py
+++ b/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/makeIC.py
@@ -1,22 +1,22 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
- #               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) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+#               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/>.
+#
+##############################################################################
 
 import h5py
 import numpy as np
@@ -24,27 +24,27 @@ import numpy as np
 # Generates a swift IC file for the Kelvin-Helmholtz test in a periodic box
 
 # Constants
-R_earth = 6371000    # Earth radius
-R_jupiter = 11.2089 * R_earth   # Jupiter radius
+R_earth = 6371000  # Earth radius
+R_jupiter = 11.2089 * R_earth  # Jupiter radius
 
 # Parameters
-N2_l     = 128       # Particles along one edge in the low-density region
-N2_depth = 18        # Particles in z direction in low-density region
-matID1 = 304         # Central region material ID: AQUA
-matID2 = 307         # Outskirts material ID: CD21 H--He
-P1    = 3.2e+12      # Central region pressure
-P2    = 3.2e+12      # Outskirts pressure
-u1    = 283591514    # Central region specific internal energy
-u2    = 804943158    # Outskirts specific internal energy
-rho1_approx =  9000  # Central region density. Readjusted later
-rho2 =  3500         # Outskirts density
-boxsize_l = R_jupiter     # size of simulation box in x and y dimension
-v1   = boxsize_l / 10000  # Central region velocity
-v2   = -boxsize_l / 10000 # Outskirts velocity
-boxsize_depth = boxsize_l * N2_depth / N2_l # size of simulation box in z dimension
+N2_l = 128  # Particles along one edge in the low-density region
+N2_depth = 18  # Particles in z direction in low-density region
+matID1 = 304  # Central region material ID: AQUA
+matID2 = 307  # Outskirts material ID: CD21 H--He
+P1 = 3.2e12  # Central region pressure
+P2 = 3.2e12  # Outskirts pressure
+u1 = 283591514  # Central region specific internal energy
+u2 = 804943158  # Outskirts specific internal energy
+rho1_approx = 9000  # Central region density. Readjusted later
+rho2 = 3500  # Outskirts density
+boxsize_l = R_jupiter  # size of simulation box in x and y dimension
+v1 = boxsize_l / 10000  # Central region velocity
+v2 = -boxsize_l / 10000  # Outskirts velocity
+boxsize_depth = boxsize_l * N2_depth / N2_l  # size of simulation box in z dimension
 mass = rho2 * (boxsize_l * boxsize_l * boxsize_depth) / (N2_l * N2_l * N2_depth)
 fileOutputName = "kelvin_helmholtz.hdf5"
-#---------------------------------------------------
+# ---------------------------------------------------
 
 # Start by calculating N1_l and rho1
 numPart2 = N2_l * N2_l * N2_depth
@@ -87,23 +87,27 @@ for i in range(N1_depth):
     for j in range(N1_l):
         for k in range(N1_l):
             index = i * N1_l**2 + j * N1_l + k
-            A2_coords1[index, 0] = (j / float(N1_l) + 1. / (2. * N1_l)) * boxsize_l
-            A2_coords1[index, 1] = (k / float(N1_l) + 1. / (2. * N1_l)) * boxsize_l
-            A2_coords1[index, 2] = (i / float(N1_depth) + 1. / (2. * N1_depth)) * boxsize_depth
+            A2_coords1[index, 0] = (j / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_l
+            A2_coords1[index, 1] = (k / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_l
+            A2_coords1[index, 2] = (
+                i / float(N1_depth) + 1.0 / (2.0 * N1_depth)
+            ) * boxsize_depth
 
 # Particles in the outskirts
 for i in range(N2_depth):
     for j in range(N2_l):
         for k in range(N2_l):
             index = i * N2_l**2 + j * N2_l + k
-            A2_coords2[index, 0] = (j / float(N2_l) + 1. / (2. * N2_l)) * boxsize_l
-            A2_coords2[index, 1] = (k / float(N2_l) + 1. / (2. * N2_l)) * boxsize_l
-            A2_coords2[index, 2] = (i / float(N2_depth) + 1. / (2. * N2_depth)) * boxsize_depth
+            A2_coords2[index, 0] = (j / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_l
+            A2_coords2[index, 1] = (k / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_l
+            A2_coords2[index, 2] = (
+                i / float(N2_depth) + 1.0 / (2.0 * N2_depth)
+            ) * boxsize_depth
 
 
 # Masks for the particles to be selected for the outer and inner regions
-mask1 = abs(A2_coords1[:,1] - 0.5 * boxsize_l) < 0.25 * boxsize_l
-mask2 = abs(A2_coords2[:,1] - 0.5 * boxsize_l) > 0.25 * boxsize_l
+mask1 = abs(A2_coords1[:, 1] - 0.5 * boxsize_l) < 0.25 * boxsize_l
+mask2 = abs(A2_coords2[:, 1] - 0.5 * boxsize_l) > 0.25 * boxsize_l
 
 # The positions of the particles are now selected
 # and the placement of the lattices are adjusted to give appropriate interfaces
@@ -116,16 +120,22 @@ pcl_separation_2 = np.cbrt(mass / rho2)
 boundary_separation = 0.5 * (pcl_separation_1 + pcl_separation_2)
 
 # Shift all the "inside" particles to get boundary_separation across the bottom interface
-min_y_inside = np.min(A2_coords_inside[:,1])
-max_y_outside_bot = np.max(A2_coords_outside[A2_coords_outside[:,1]-0.5 * boxsize_l < -0.25 * boxsize_l, 1])
+min_y_inside = np.min(A2_coords_inside[:, 1])
+max_y_outside_bot = np.max(
+    A2_coords_outside[A2_coords_outside[:, 1] - 0.5 * boxsize_l < -0.25 * boxsize_l, 1]
+)
 shift_distance_bot = boundary_separation - (min_y_inside - max_y_outside_bot)
-A2_coords_inside[:,1] += shift_distance_bot
+A2_coords_inside[:, 1] += shift_distance_bot
 
 # Shift the top section of the "outside" particles to get boundary_separation across the top interface
-max_y_inside = np.max(A2_coords_inside[:,1])
-min_y_outside_top = np.min(A2_coords_outside[A2_coords_outside[:,1]-0.5 * boxsize_l > 0.25 * boxsize_l, 1])
+max_y_inside = np.max(A2_coords_inside[:, 1])
+min_y_outside_top = np.min(
+    A2_coords_outside[A2_coords_outside[:, 1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1]
+)
 shift_distance_top = boundary_separation - (min_y_outside_top - max_y_inside)
-A2_coords_outside[A2_coords_outside[:,1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1] += shift_distance_top
+A2_coords_outside[
+    A2_coords_outside[:, 1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1
+] += shift_distance_top
 
 # Adjust box size in y direction based on the shifting of the lattices.
 new_box_y = boxsize_l + shift_distance_top
@@ -143,14 +153,16 @@ A1_ids = np.linspace(1, numPart, numPart)
 
 # Finally add the velocity perturbation
 vel_perturb_factor = 0.01 * (v1 - v2)
-A2_vel[:,1] = vel_perturb_factor * np.sin(2 * np.pi * A2_coords[:,0] / (0.5 * boxsize_l))
+A2_vel[:, 1] = vel_perturb_factor * np.sin(
+    2 * np.pi * A2_coords[:, 0] / (0.5 * boxsize_l)
+)
 
 # Write ICs to file
 with h5py.File(fileOutputName, "w") as f:
     # Header
     grp = f.create_group("/Header")
     grp.attrs["BoxSize"] = [boxsize_l, new_box_y, boxsize_depth]
-    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
@@ -159,29 +171,29 @@ with h5py.File(fileOutputName, "w") as f:
     grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
     grp.attrs["Dimension"] = 3
 
-    #Units
+    # Units
     grp = f.create_group("/Units")
-    grp.attrs["Unit length in cgs (U_L)"] = 100.
-    grp.attrs["Unit mass in cgs (U_M)"] = 1000.
-    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)"] = 100.0
+    grp.attrs["Unit mass in cgs (U_M)"] = 1000.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 = f.create_group("/PartType0")
-    ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+    ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
     ds[()] = A2_coords
-    ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+    ds = grp.create_dataset("Velocities", (numPart, 3), "f")
     ds[()] = A2_vel
-    ds = grp.create_dataset('Masses', (numPart, 1), 'f')
-    ds[()] = A1_m.reshape((numPart,1))
-    ds = grp.create_dataset('Density', (numPart,1), 'f')
-    ds[()] = A1_rho.reshape((numPart,1))
-    ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-    ds[()] = A1_h.reshape((numPart,1))
-    ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-    ds[()] = A1_u.reshape((numPart,1))
-    ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
-    ds[()] = A1_ids.reshape((numPart,1))
-    ds = grp.create_dataset('MaterialIDs', (numPart,1), 'i')
-    ds[()] = A1_mat.reshape((numPart,1))
+    ds = grp.create_dataset("Masses", (numPart, 1), "f")
+    ds[()] = A1_m.reshape((numPart, 1))
+    ds = grp.create_dataset("Density", (numPart, 1), "f")
+    ds[()] = A1_rho.reshape((numPart, 1))
+    ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+    ds[()] = A1_h.reshape((numPart, 1))
+    ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+    ds[()] = A1_u.reshape((numPart, 1))
+    ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+    ds[()] = A1_ids.reshape((numPart, 1))
+    ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
+    ds[()] = A1_mat.reshape((numPart, 1))
diff --git a/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/plotModeGrowth.py b/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/plotModeGrowth.py
index 5ff9211bd4031c14463c75d7ce237f6b304ee5bd..3785996dc39a1f31d92e1ed9351b5ff20b06f875 100644
--- a/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/plotModeGrowth.py
+++ b/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/plotModeGrowth.py
@@ -28,6 +28,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def calculate_mode_growth(snap, vy_init_amp, wavelength):
 
     # Load data
@@ -40,46 +41,72 @@ def calculate_mode_growth(snap, vy_init_amp, wavelength):
         A1_h = f["/PartType0/SmoothingLengths"][:] / boxsize_l
 
     # Adjust positions
-    A2_pos[:,0] += 0.5
-    A2_pos[:,1] += 0.5
+    A2_pos[:, 0] += 0.5
+    A2_pos[:, 1] += 0.5
 
     # Masks to select the upper and lower halfs of the simulations
-    mask_up = A2_pos[:,1] >= 0.5
-    mask_down = A2_pos[:,1] < 0.5
+    mask_up = A2_pos[:, 1] >= 0.5
+    mask_down = A2_pos[:, 1] < 0.5
 
     # McNally et al. 2012 Eqn. 10
     s = np.empty(len(A1_h))
-    s[mask_down] = A2_vel[mask_down, 1] * A1_h[mask_down]**3 * np.sin(2 * np.pi * A2_pos[mask_down, 0] / wavelength) * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    s[mask_up] = A2_vel[mask_up, 1] * A1_h[mask_up]**3 * np.sin(2 * np.pi * A2_pos[mask_up, 0] / wavelength) * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    s[mask_down] = (
+        A2_vel[mask_down, 1]
+        * A1_h[mask_down] ** 3
+        * np.sin(2 * np.pi * A2_pos[mask_down, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
+    )
+    s[mask_up] = (
+        A2_vel[mask_up, 1]
+        * A1_h[mask_up] ** 3
+        * np.sin(2 * np.pi * A2_pos[mask_up, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    )
 
     # McNally et al. 2012 Eqn. 11
     c = np.empty(len(A1_h))
-    c[mask_down] = A2_vel[mask_down, 1] * A1_h[mask_down]**3 * np.cos(2 * np.pi * A2_pos[mask_down, 0] / wavelength) * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    c[mask_up] = A2_vel[mask_up, 1] * A1_h[mask_up]**3 * np.cos(2 * np.pi * A2_pos[mask_up, 0] / wavelength) * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    c[mask_down] = (
+        A2_vel[mask_down, 1]
+        * A1_h[mask_down] ** 3
+        * np.cos(2 * np.pi * A2_pos[mask_down, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
+    )
+    c[mask_up] = (
+        A2_vel[mask_up, 1]
+        * A1_h[mask_up] ** 3
+        * np.cos(2 * np.pi * A2_pos[mask_up, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    )
 
     # McNally et al. 2012 Eqn. 12
     d = np.empty(len(A1_h))
-    d[mask_down] = A1_h[mask_down]**3 * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
-    d[mask_up] = A1_h[mask_up]**3 * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    d[mask_down] = A1_h[mask_down] ** 3 * np.exp(
+        -2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength
+    )
+    d[mask_up] = A1_h[mask_up] ** 3 * np.exp(
+        -2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength
+    )
 
     # McNally et al. 2012 Eqn. 13
-    M = 2 * np.sqrt((np.sum(s)/ np.sum(d))**2 + (np.sum(c)/ np.sum(d))**2)
+    M = 2 * np.sqrt((np.sum(s) / np.sum(d)) ** 2 + (np.sum(c) / np.sum(d)) ** 2)
 
     return M
 
 
-
 if __name__ == "__main__":
 
     # Simulation paramerters for nomralisation of mode growth
-    vy_init_amp = 2e-6 # Initial amplitude of y velocity perturbation in units of the boxsize per second
-    wavelength = 0.5   # wavelength of initial perturbation in units of the boxsize
+    vy_init_amp = 2e-6  # Initial amplitude of y velocity perturbation in units of the boxsize per second
+    wavelength = 0.5  # wavelength of initial perturbation in units of the boxsize
 
     # Use Gridspec to set up figure
     n_gs_ax = 40
     ax_len = 5
     fig = plt.figure(figsize=(9, 6))
-    gs = mpl.gridspec.GridSpec(nrows=40, ncols=60,)
+    gs = mpl.gridspec.GridSpec(
+        nrows=40,
+        ncols=60,
+    )
     ax = plt.subplot(gs[:, :])
 
     # Snapshots and corresponding times
@@ -98,8 +125,8 @@ if __name__ == "__main__":
     ax.set_ylabel(r"$\log( \, M \; / \; M^{}_0 \, )$", fontsize=18)
     ax.set_xlabel(r"$t \; / \; \tau^{}_{\rm KH}$", fontsize=18)
     ax.minorticks_on()
-    ax.tick_params(which='major', direction="in")
-    ax.tick_params(which='minor', direction="in")
+    ax.tick_params(which="major", direction="in")
+    ax.tick_params(which="minor", direction="in")
     ax.tick_params(labelsize=14)
 
     plt.savefig("mode_growth.png", dpi=300, bbox_inches="tight")
diff --git a/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/plotSnapshots.py b/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/plotSnapshots.py
index 18df751c1040b061ecc00c41afd2847de935a68d..818be1849ceb97d4461915b0159627b73c399545 100644
--- a/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/plotSnapshots.py
+++ b/examples/Planetary/KelvinHelmholtz_JupiterLike_3D/plotSnapshots.py
@@ -27,6 +27,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def make_axes():
     # Use Gridspec to set up figure
     n_gs_ax = 40
@@ -48,11 +49,33 @@ def make_axes():
     )
 
     ax_0 = plt.subplot(gs[:n_gs_ax, :n_gs_ax])
-    ax_1 = plt.subplot(gs[:n_gs_ax, n_gs_ax+n_gs_ax_gap:2*n_gs_ax+n_gs_ax_gap])
-    ax_2 = plt.subplot(gs[:n_gs_ax, 2*n_gs_ax+2*n_gs_ax_gap:3*n_gs_ax+2*n_gs_ax_gap])
+    ax_1 = plt.subplot(gs[:n_gs_ax, n_gs_ax + n_gs_ax_gap : 2 * n_gs_ax + n_gs_ax_gap])
+    ax_2 = plt.subplot(
+        gs[:n_gs_ax, 2 * n_gs_ax + 2 * n_gs_ax_gap : 3 * n_gs_ax + 2 * n_gs_ax_gap]
+    )
 
-    cax_0 = plt.subplot(gs[:int(0.5*n_gs_ax)-1, 3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap:3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap+n_gs_cbar])
-    cax_1 = plt.subplot(gs[int(0.5*n_gs_ax)+1:, 3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap:3*n_gs_ax+2*n_gs_ax_gap+n_gs_cbar_gap+n_gs_cbar])
+    cax_0 = plt.subplot(
+        gs[
+            : int(0.5 * n_gs_ax) - 1,
+            3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
+    cax_1 = plt.subplot(
+        gs[
+            int(0.5 * n_gs_ax) + 1 :,
+            3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
 
     axs = [ax_0, ax_1, ax_2]
     caxs = [cax_0, cax_1]
@@ -67,8 +90,8 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
 
     with h5py.File(snap_file, "r") as f:
         # Units from file metadata to SI
-        m=float(f["Units"].attrs["Unit mass in cgs (U_M)"][0]) * 1e-3
-        l=float(f["Units"].attrs["Unit length in cgs (U_L)"][0]) * 1e-2
+        m = float(f["Units"].attrs["Unit mass in cgs (U_M)"][0]) * 1e-3
+        l = float(f["Units"].attrs["Unit length in cgs (U_L)"][0]) * 1e-2
 
         boxsize_l = f["Header"].attrs["BoxSize"][0] * l
         A1_x = f["/PartType0/Coordinates"][:, 0] * l
@@ -137,21 +160,20 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
     ax.set_ylim((0, boxsize_l))
 
 
-
 if __name__ == "__main__":
 
     # Set colormap
-    cmap1 = plt.get_cmap('winter_r')
+    cmap1 = plt.get_cmap("winter_r")
     mat_id1 = 304
     rho_min1 = 6800
     rho_max1 = 8700
-    norm1 = mpl.colors.Normalize(vmin=rho_min1,vmax=rho_max1)
+    norm1 = mpl.colors.Normalize(vmin=rho_min1, vmax=rho_max1)
 
-    cmap2 = plt.get_cmap('autumn')
+    cmap2 = plt.get_cmap("autumn")
     mat_id2 = 307
     rho_min2 = 3450
     rho_max2 = 4050
-    norm2 = mpl.colors.Normalize(vmin=rho_min2,vmax=rho_max2)
+    norm2 = mpl.colors.Normalize(vmin=rho_min2, vmax=rho_max2)
 
     # Generate axes
     axs, caxs = make_axes()
@@ -166,7 +188,14 @@ if __name__ == "__main__":
         time = times[i]
 
         plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2)
-        ax.text(0.5,-0.1,r"$t =\;$"+time+r"$\, \tau^{}_{\rm KH}$", horizontalalignment='center',size=18, transform=ax.transAxes)
+        ax.text(
+            0.5,
+            -0.1,
+            r"$t =\;$" + time + r"$\, \tau^{}_{\rm KH}$",
+            horizontalalignment="center",
+            size=18,
+            transform=ax.transAxes,
+        )
 
     # Colour bar
     sm1 = plt.cm.ScalarMappable(cmap=cmap1, norm=norm1)
diff --git a/examples/Planetary/RayleighTaylor_EarthLike_3D/makeIC.py b/examples/Planetary/RayleighTaylor_EarthLike_3D/makeIC.py
index da4d7ee11c9e7e22752ed2c0a84e99fac4a84efe..07a4c16b82afe2a5a756b80d37b59925c77c5958 100644
--- a/examples/Planetary/RayleighTaylor_EarthLike_3D/makeIC.py
+++ b/examples/Planetary/RayleighTaylor_EarthLike_3D/makeIC.py
@@ -1,58 +1,67 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
- #               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) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+#               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/>.
+#
+##############################################################################
 
 import h5py
 import numpy as np
 import woma
 
 # Load EoS tables
-woma.load_eos_tables(["ANEOS_forsterite","ANEOS_Fe85Si15"])
+woma.load_eos_tables(["ANEOS_forsterite", "ANEOS_Fe85Si15"])
 
 # Generates a swift IC file for the Rayleigh-Taylor instability test
 
 # Constants
-R_earth = 6371000 # Earth radius
+R_earth = 6371000  # Earth radius
 
 # Parameters
-N2_l = 64             # Number of particles along one edge in lower region
-N2_depth = 18         # Number of particles along in z dimension in lower region
-matID1 = 402          # Upper region material ID: ANEOS Fe85Si15
-matID2 = 400          # Lower region material ID: ANEOS forsterite
-rho1_approx = 10000   # Approximate density of upper region. To be recalculated
-rho2 = 5000           # Density of lower region
-g = -9.91             # Constant gravitational acceleration
-P0 = 1.2e+11          # Pressure at interface
+N2_l = 64  # Number of particles along one edge in lower region
+N2_depth = 18  # Number of particles along in z dimension in lower region
+matID1 = 402  # Upper region material ID: ANEOS Fe85Si15
+matID2 = 400  # Lower region material ID: ANEOS forsterite
+rho1_approx = 10000  # Approximate density of upper region. To be recalculated
+rho2 = 5000  # Density of lower region
+g = -9.91  # Constant gravitational acceleration
+P0 = 1.2e11  # Pressure at interface
 boxsize_factor = 0.1 * R_earth
 dv = 0.00025 * boxsize_factor  # Size of velocity perturbation
-boxsize_xy = [0.5 * boxsize_factor, 1.0 * boxsize_factor]  # Size of the box in x and y dimensions
-boxsize_depth = boxsize_xy[0] * N2_depth / N2_l # Size of simulation box in z dimension
-fixed_region = [0.05 * boxsize_factor, 0.95 * boxsize_factor]       # y-range of non fixed_region particles
-perturbation_region = [0.3 * boxsize_factor, 0.7 * boxsize_factor]  # y-range for the velocity perturbation
+boxsize_xy = [
+    0.5 * boxsize_factor,
+    1.0 * boxsize_factor,
+]  # Size of the box in x and y dimensions
+boxsize_depth = boxsize_xy[0] * N2_depth / N2_l  # Size of simulation box in z dimension
+fixed_region = [
+    0.05 * boxsize_factor,
+    0.95 * boxsize_factor,
+]  # y-range of non fixed_region particles
+perturbation_region = [
+    0.3 * boxsize_factor,
+    0.7 * boxsize_factor,
+]  # y-range for the velocity perturbation
 fileOutputName = "rayleigh_taylor.hdf5"
-#---------------------------------------------------
+# ---------------------------------------------------
 
 # Start by generating grids of particles of the two densities
 numPart2 = N2_l * N2_l * N2_depth
 numPart1 = int(numPart2 / rho2 * rho1_approx)
 N1_l = int(np.cbrt(boxsize_xy[0] * numPart1 / boxsize_depth))
-N1_l -= N1_l % 4 # Make RT symmetric across centre of both instability regions
+N1_l -= N1_l % 4  # Make RT symmetric across centre of both instability regions
 N1_depth = int(boxsize_depth * N1_l / boxsize_xy[0])
 numPart1 = int(N1_l * N1_l * N1_depth)
 numPart = numPart2 + numPart1
@@ -89,13 +98,22 @@ for i in range(N1_depth):
     for j in range(N1_l):
         for k in range(N1_l):
             index = i * N1_l**2 + j * N1_l + k
-            A2_coords1[index, 0] = (j / float(N1_l) + 1. / (2. * N1_l)) * boxsize_xy[0]
-            A2_coords1[index, 1] = (k / float(N1_l) + 1. / (2. * N1_l)) * (0.5 * boxsize_xy[1]) + 0.5 * boxsize_xy[1]
-            A2_coords1[index, 2] = (i / float(N1_depth) + 1. / (2. * N1_depth)) * boxsize_depth
+            A2_coords1[index, 0] = (j / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_xy[
+                0
+            ]
+            A2_coords1[index, 1] = (k / float(N1_l) + 1.0 / (2.0 * N1_l)) * (
+                0.5 * boxsize_xy[1]
+            ) + 0.5 * boxsize_xy[1]
+            A2_coords1[index, 2] = (
+                i / float(N1_depth) + 1.0 / (2.0 * N1_depth)
+            ) * boxsize_depth
             A1_rho1[index] = rho1
 
             # If in top and bottom where particles are fixed
-            if A2_coords1[index, 1] < fixed_region[0] or A2_coords1[index, 1] > fixed_region[1]:
+            if (
+                A2_coords1[index, 1] < fixed_region[0]
+                or A2_coords1[index, 1] > fixed_region[1]
+            ):
                 A1_ids[index] = boundary_particles
                 boundary_particles += 1
 
@@ -105,22 +123,34 @@ for i in range(N2_depth):
     for j in range(N2_l):
         for k in range(N2_l):
             index = i * N2_l**2 + j * N2_l + k
-            A2_coords2[index, 0] = (j / float(N2_l) + 1. / (2. * N2_l)) * boxsize_xy[0]
-            A2_coords2[index, 1] = (k / float(N2_l) + 1. / (2. * N2_l)) * (0.5 * boxsize_xy[1])
-            A2_coords2[index, 2] = (i / float(N2_depth) + 1. / (2. * N2_depth)) * boxsize_depth
+            A2_coords2[index, 0] = (j / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_xy[
+                0
+            ]
+            A2_coords2[index, 1] = (k / float(N2_l) + 1.0 / (2.0 * N2_l)) * (
+                0.5 * boxsize_xy[1]
+            )
+            A2_coords2[index, 2] = (
+                i / float(N2_depth) + 1.0 / (2.0 * N2_depth)
+            ) * boxsize_depth
             A1_rho2[index] = rho2
 
             # If in top and bottom where particles are fixed
-            if A2_coords2[index, 1] < fixed_region[0] or A2_coords2[index, 1] > fixed_region[1]:
+            if (
+                A2_coords2[index, 1] < fixed_region[0]
+                or A2_coords2[index, 1] > fixed_region[1]
+            ):
                 A1_ids[index + numPart1] = boundary_particles
                 boundary_particles += 1
 
 print(
-    "You need to compile the code with " "--enable-boundary-particles=%i" % boundary_particles
+    "You need to compile the code with "
+    "--enable-boundary-particles=%i" % boundary_particles
 )
 
 # Set IDs of non-boundary particles
-A1_ids[A1_ids == 0] = np.linspace(boundary_particles, numPart, numPart - boundary_particles + 1)
+A1_ids[A1_ids == 0] = np.linspace(
+    boundary_particles, numPart, numPart - boundary_particles + 1
+)
 
 # The placement of the lattices are now adjusted to give appropriate interfaces
 # Calculate the separation of particles across the density discontinuity
@@ -129,16 +159,16 @@ pcl_separation_1 = np.cbrt(mass / rho1)
 boundary_separation = 0.5 * (pcl_separation_2 + pcl_separation_1)
 
 # Shift top lattice
-min_y1 = min(A2_coords1[:,1])
+min_y1 = min(A2_coords1[:, 1])
 max_y2 = max(A2_coords2[:, 1])
 shift_distance = boundary_separation - (min_y1 - max_y2)
-A2_coords1[:,1] += shift_distance
+A2_coords1[:, 1] += shift_distance
 
 # Calculate internal energies
 A1_P1 = P0 + g * A1_rho1 * (A2_coords1[:, 1] - 0.5 * boxsize_xy[1])
 A1_P2 = P0 + g * A1_rho2 * (A2_coords2[:, 1] - 0.5 * boxsize_xy[1])
-A1_u1 =  woma.A1_Z_rho_Y(A1_rho1, A1_P1, A1_mat1, Z_choice="u", Y_choice="P")
-A1_u2 =  woma.A1_Z_rho_Y(A1_rho2, A1_P2, A1_mat2, Z_choice="u", Y_choice="P")
+A1_u1 = woma.A1_Z_rho_Y(A1_rho1, A1_P1, A1_mat1, Z_choice="u", Y_choice="P")
+A1_u2 = woma.A1_Z_rho_Y(A1_rho2, A1_P2, A1_mat2, Z_choice="u", Y_choice="P")
 
 # Now the two lattices can be combined
 A2_coords = np.append(A2_coords1, A2_coords2, axis=0)
@@ -151,18 +181,26 @@ A1_h = np.append(A1_h1, A1_h2, axis=0)
 
 # Add velocity perturbation
 mask_perturb = np.logical_and(
-                               A2_coords[:, 1] > perturbation_region[0],
-                               A2_coords[:, 1] < perturbation_region[1],
-                              )
-A2_vel[mask_perturb, 1] = dv * (1 + np.cos(8 * np.pi * (A2_coords[mask_perturb, 0] / (boxsize_factor) + 0.25))) * (1 + np.cos(5 * np.pi * (A2_coords[mask_perturb, 1] / (boxsize_factor) - 0.5)))
+    A2_coords[:, 1] > perturbation_region[0],
+    A2_coords[:, 1] < perturbation_region[1],
+)
+A2_vel[mask_perturb, 1] = (
+    dv
+    * (1 + np.cos(8 * np.pi * (A2_coords[mask_perturb, 0] / (boxsize_factor) + 0.25)))
+    * (1 + np.cos(5 * np.pi * (A2_coords[mask_perturb, 1] / (boxsize_factor) - 0.5)))
+)
 
 
 # Write ICs to file
 with h5py.File(fileOutputName, "w") as f:
     # Header
     grp = f.create_group("/Header")
-    grp.attrs["BoxSize"] = [boxsize_xy[0], boxsize_xy[1] + shift_distance, boxsize_depth]
-    grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+    grp.attrs["BoxSize"] = [
+        boxsize_xy[0],
+        boxsize_xy[1] + shift_distance,
+        boxsize_depth,
+    ]
+    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
@@ -171,29 +209,29 @@ with h5py.File(fileOutputName, "w") as f:
     grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
     grp.attrs["Dimension"] = 3
 
-    #Units
+    # Units
     grp = f.create_group("/Units")
-    grp.attrs["Unit length in cgs (U_L)"] = 100.
-    grp.attrs["Unit mass in cgs (U_M)"] = 1000.
-    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)"] = 100.0
+    grp.attrs["Unit mass in cgs (U_M)"] = 1000.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 = f.create_group("/PartType0")
-    ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+    ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
     ds[()] = A2_coords
-    ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+    ds = grp.create_dataset("Velocities", (numPart, 3), "f")
     ds[()] = A2_vel
-    ds = grp.create_dataset('Masses', (numPart, 1), 'f')
-    ds[()] = A1_m.reshape((numPart,1))
-    ds = grp.create_dataset('Density', (numPart,1), 'f')
-    ds[()] = A1_rho.reshape((numPart,1))
-    ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-    ds[()] = A1_h.reshape((numPart,1))
-    ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-    ds[()] = A1_u.reshape((numPart,1))
-    ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
-    ds[()] = A1_ids.reshape((numPart,1))
-    ds = grp.create_dataset('MaterialIDs', (numPart,1), 'i')
-    ds[()] = A1_mat.reshape((numPart,1))
+    ds = grp.create_dataset("Masses", (numPart, 1), "f")
+    ds[()] = A1_m.reshape((numPart, 1))
+    ds = grp.create_dataset("Density", (numPart, 1), "f")
+    ds[()] = A1_rho.reshape((numPart, 1))
+    ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+    ds[()] = A1_h.reshape((numPart, 1))
+    ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+    ds[()] = A1_u.reshape((numPart, 1))
+    ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+    ds[()] = A1_ids.reshape((numPart, 1))
+    ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
+    ds[()] = A1_mat.reshape((numPart, 1))
diff --git a/examples/Planetary/RayleighTaylor_EarthLike_3D/plotSnapshots.py b/examples/Planetary/RayleighTaylor_EarthLike_3D/plotSnapshots.py
index edb29f98e1892a759dc4d7cb16755c171122b435..fee140e5fb5803de9eaadf787e59b4db95aaf751 100644
--- a/examples/Planetary/RayleighTaylor_EarthLike_3D/plotSnapshots.py
+++ b/examples/Planetary/RayleighTaylor_EarthLike_3D/plotSnapshots.py
@@ -27,6 +27,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def make_axes():
     # Use Gridspec to set up figure
     n_gs_ax_x = 40
@@ -50,11 +51,38 @@ def make_axes():
     )
 
     ax_0 = plt.subplot(gs[:n_gs_ax_y, :n_gs_ax_x])
-    ax_1 = plt.subplot(gs[:n_gs_ax_y, n_gs_ax_x+n_gs_ax_gap:2*n_gs_ax_x+n_gs_ax_gap])
-    ax_2 = plt.subplot(gs[:n_gs_ax_y, 2*n_gs_ax_x+2*n_gs_ax_gap:3*n_gs_ax_x+2*n_gs_ax_gap])
+    ax_1 = plt.subplot(
+        gs[:n_gs_ax_y, n_gs_ax_x + n_gs_ax_gap : 2 * n_gs_ax_x + n_gs_ax_gap]
+    )
+    ax_2 = plt.subplot(
+        gs[
+            :n_gs_ax_y,
+            2 * n_gs_ax_x + 2 * n_gs_ax_gap : 3 * n_gs_ax_x + 2 * n_gs_ax_gap,
+        ]
+    )
 
-    cax_0 = plt.subplot(gs[:int(0.5*n_gs_ax_y)-1, 3*n_gs_ax_x+2*n_gs_ax_gap+n_gs_cbar_gap:3*n_gs_ax_x+2*n_gs_ax_gap+n_gs_cbar_gap+n_gs_cbar])
-    cax_1 = plt.subplot(gs[int(0.5*n_gs_ax_y)+1:, 3*n_gs_ax_x+2*n_gs_ax_gap+n_gs_cbar_gap:3*n_gs_ax_x+2*n_gs_ax_gap+n_gs_cbar_gap+n_gs_cbar])
+    cax_0 = plt.subplot(
+        gs[
+            : int(0.5 * n_gs_ax_y) - 1,
+            3 * n_gs_ax_x
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax_x
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
+    cax_1 = plt.subplot(
+        gs[
+            int(0.5 * n_gs_ax_y) + 1 :,
+            3 * n_gs_ax_x
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax_x
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
 
     axs = [ax_0, ax_1, ax_2]
     caxs = [cax_0, cax_1]
@@ -69,8 +97,8 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
 
     with h5py.File(snap_file, "r") as f:
         # Units from file metadata to SI
-        m=float(f["Units"].attrs["Unit mass in cgs (U_M)"][0]) * 1e-3
-        l=float(f["Units"].attrs["Unit length in cgs (U_L)"][0]) * 1e-2
+        m = float(f["Units"].attrs["Unit mass in cgs (U_M)"][0]) * 1e-3
+        l = float(f["Units"].attrs["Unit length in cgs (U_L)"][0]) * 1e-2
 
         boxsize_x = f["Header"].attrs["BoxSize"][0] * l
         boxsize_y = f["Header"].attrs["BoxSize"][1] * l
@@ -140,21 +168,20 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
     ax.set_ylim((0.05 * boxsize_y, 0.95 * boxsize_y))
 
 
-
 if __name__ == "__main__":
 
     # Set colormap
-    cmap1 = plt.get_cmap('YlOrRd')
+    cmap1 = plt.get_cmap("YlOrRd")
     mat_id1 = 402
     rho_min1 = 7950
     rho_max1 = 10050
-    norm1 = mpl.colors.Normalize(vmin=rho_min1,vmax=rho_max1)
+    norm1 = mpl.colors.Normalize(vmin=rho_min1, vmax=rho_max1)
 
-    cmap2 = plt.get_cmap('Blues_r')
+    cmap2 = plt.get_cmap("Blues_r")
     mat_id2 = 400
     rho_min2 = 4950
     rho_max2 = 5550
-    norm2 = mpl.colors.Normalize(vmin=rho_min2,vmax=rho_max2)
+    norm2 = mpl.colors.Normalize(vmin=rho_min2, vmax=rho_max2)
 
     # Generate axes
     axs, caxs = make_axes()
@@ -169,7 +196,14 @@ if __name__ == "__main__":
         time = times[i]
 
         plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2)
-        ax.text(0.5,-0.05,r"$t =\;$"+time+r"$\,$s", horizontalalignment='center',size=18, transform=ax.transAxes)
+        ax.text(
+            0.5,
+            -0.05,
+            r"$t =\;$" + time + r"$\,$s",
+            horizontalalignment="center",
+            size=18,
+            transform=ax.transAxes,
+        )
 
     # Colour bar
     sm1 = plt.cm.ScalarMappable(cmap=cmap1, norm=norm1)
diff --git a/examples/Planetary/RayleighTaylor_IdealGas_3D/makeIC.py b/examples/Planetary/RayleighTaylor_IdealGas_3D/makeIC.py
index 294b00f2d9ee086d953c214144a12e924ab6ce90..1ac7d821b1bc99107124be1999de3065ee095499 100644
--- a/examples/Planetary/RayleighTaylor_IdealGas_3D/makeIC.py
+++ b/examples/Planetary/RayleighTaylor_IdealGas_3D/makeIC.py
@@ -1,22 +1,22 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
- #               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) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+#               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/>.
+#
+##############################################################################
 
 import h5py
 import numpy as np
@@ -24,28 +24,37 @@ import numpy as np
 # Generates a swift IC file for the Rayleigh-Taylor instability test
 
 # Parameters
-N2_l = 64         # Number of particles along one edge in lower region
-N2_depth = 18     # Number of particles along in z dimension in lower region
-gamma = 7.0 / 5.0 # Gas adiabatic index
-matID1 = 0        # Upper region material ID: ideal gas
-matID2 = 0        # Lower region material ID: ideal gas
-rho1_approx = 2   # Approximate density of upper region. To be recalculated
-rho2 = 1          # Density of lower region
-g = -0.5          # Constant gravitational acceleration
-dv = 0.025        # Size of velocity perturbation
+N2_l = 64  # Number of particles along one edge in lower region
+N2_depth = 18  # Number of particles along in z dimension in lower region
+gamma = 7.0 / 5.0  # Gas adiabatic index
+matID1 = 0  # Upper region material ID: ideal gas
+matID2 = 0  # Lower region material ID: ideal gas
+rho1_approx = 2  # Approximate density of upper region. To be recalculated
+rho2 = 1  # Density of lower region
+g = -0.5  # Constant gravitational acceleration
+dv = 0.025  # Size of velocity perturbation
 boxsize_factor = 1
-boxsize_xy = [0.5 * boxsize_factor, 1.0 * boxsize_factor]  # Size of the box in x and y dimensions
-boxsize_depth = boxsize_xy[0] * N2_depth / N2_l # Size of simulation box in z dimension
-fixed_region = [0.05 * boxsize_factor, 0.95 * boxsize_factor]       # y-range of non fixed_region particles
-perturbation_region = [0.3 * boxsize_factor, 0.7 * boxsize_factor]  # y-range for the velocity perturbation
+boxsize_xy = [
+    0.5 * boxsize_factor,
+    1.0 * boxsize_factor,
+]  # Size of the box in x and y dimensions
+boxsize_depth = boxsize_xy[0] * N2_depth / N2_l  # Size of simulation box in z dimension
+fixed_region = [
+    0.05 * boxsize_factor,
+    0.95 * boxsize_factor,
+]  # y-range of non fixed_region particles
+perturbation_region = [
+    0.3 * boxsize_factor,
+    0.7 * boxsize_factor,
+]  # y-range for the velocity perturbation
 fileOutputName = "rayleigh_taylor.hdf5"
-#---------------------------------------------------
+# ---------------------------------------------------
 
 # Start by generating grids of particles of the two densities
 numPart2 = N2_l * N2_l * N2_depth
 numPart1 = int(numPart2 / rho2 * rho1_approx)
 N1_l = int(np.cbrt(boxsize_xy[0] * numPart1 / boxsize_depth))
-N1_l -= N1_l % 4 # Make RT symmetric across centre of both instability regions
+N1_l -= N1_l % 4  # Make RT symmetric across centre of both instability regions
 N1_depth = int(boxsize_depth * N1_l / boxsize_xy[0])
 numPart1 = int(N1_l * N1_l * N1_depth)
 numPart = numPart2 + numPart1
@@ -85,13 +94,22 @@ for i in range(N1_depth):
     for j in range(N1_l):
         for k in range(N1_l):
             index = i * N1_l**2 + j * N1_l + k
-            A2_coords1[index, 0] = (j / float(N1_l) + 1. / (2. * N1_l)) * boxsize_xy[0]
-            A2_coords1[index, 1] = (k / float(N1_l) + 1. / (2. * N1_l)) * (0.5 * boxsize_xy[1]) + 0.5 * boxsize_xy[1]
-            A2_coords1[index, 2] = (i / float(N1_depth) + 1. / (2. * N1_depth)) * boxsize_depth
+            A2_coords1[index, 0] = (j / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_xy[
+                0
+            ]
+            A2_coords1[index, 1] = (k / float(N1_l) + 1.0 / (2.0 * N1_l)) * (
+                0.5 * boxsize_xy[1]
+            ) + 0.5 * boxsize_xy[1]
+            A2_coords1[index, 2] = (
+                i / float(N1_depth) + 1.0 / (2.0 * N1_depth)
+            ) * boxsize_depth
             A1_rho1[index] = rho1
 
             # If in top and bottom where particles are fixed
-            if A2_coords1[index, 1] < fixed_region[0] or A2_coords1[index, 1] > fixed_region[1]:
+            if (
+                A2_coords1[index, 1] < fixed_region[0]
+                or A2_coords1[index, 1] > fixed_region[1]
+            ):
                 A1_ids[index] = boundary_particles
                 boundary_particles += 1
 
@@ -101,22 +119,34 @@ for i in range(N2_depth):
     for j in range(N2_l):
         for k in range(N2_l):
             index = i * N2_l**2 + j * N2_l + k
-            A2_coords2[index, 0] = (j / float(N2_l) + 1. / (2. * N2_l)) * boxsize_xy[0]
-            A2_coords2[index, 1] = (k / float(N2_l) + 1. / (2. * N2_l)) * (0.5 * boxsize_xy[1])
-            A2_coords2[index, 2] = (i / float(N2_depth) + 1. / (2. * N2_depth)) * boxsize_depth
+            A2_coords2[index, 0] = (j / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_xy[
+                0
+            ]
+            A2_coords2[index, 1] = (k / float(N2_l) + 1.0 / (2.0 * N2_l)) * (
+                0.5 * boxsize_xy[1]
+            )
+            A2_coords2[index, 2] = (
+                i / float(N2_depth) + 1.0 / (2.0 * N2_depth)
+            ) * boxsize_depth
             A1_rho2[index] = rho2
 
             # If in top and bottom where particles are fixed
-            if A2_coords2[index, 1] < fixed_region[0] or A2_coords2[index, 1] > fixed_region[1]:
+            if (
+                A2_coords2[index, 1] < fixed_region[0]
+                or A2_coords2[index, 1] > fixed_region[1]
+            ):
                 A1_ids[index + numPart1] = boundary_particles
                 boundary_particles += 1
 
 print(
-    "You need to compile the code with " "--enable-boundary-particles=%i" % boundary_particles
+    "You need to compile the code with "
+    "--enable-boundary-particles=%i" % boundary_particles
 )
 
 # Set IDs of non-boundary particles
-A1_ids[A1_ids == 0] = np.linspace(boundary_particles, numPart, numPart - boundary_particles + 1)
+A1_ids[A1_ids == 0] = np.linspace(
+    boundary_particles, numPart, numPart - boundary_particles + 1
+)
 
 # The placement of the lattices are now adjusted to give appropriate interfaces
 # Calculate the separation of particles across the density discontinuity
@@ -125,16 +155,16 @@ pcl_separation_1 = np.cbrt(mass / rho1)
 boundary_separation = 0.5 * (pcl_separation_2 + pcl_separation_1)
 
 # Shift top lattice
-min_y1 = min(A2_coords1[:,1])
+min_y1 = min(A2_coords1[:, 1])
 max_y2 = max(A2_coords2[:, 1])
 shift_distance = boundary_separation - (min_y1 - max_y2)
-A2_coords1[:,1] += shift_distance
+A2_coords1[:, 1] += shift_distance
 
 # Calculate internal energies
 A1_P1 = P0 + g * A1_rho1 * (A2_coords1[:, 1] - 0.5 * boxsize_xy[1])
 A1_P2 = P0 + g * A1_rho2 * (A2_coords2[:, 1] - 0.5 * boxsize_xy[1])
-A1_u1 =  A1_P1 / (A1_rho1 * (gamma - 1.0))
-A1_u2 =  A1_P2 / (A1_rho2 * (gamma - 1.0))
+A1_u1 = A1_P1 / (A1_rho1 * (gamma - 1.0))
+A1_u2 = A1_P2 / (A1_rho2 * (gamma - 1.0))
 
 # Now the two lattices can be combined
 A2_coords = np.append(A2_coords1, A2_coords2, axis=0)
@@ -147,17 +177,25 @@ A1_h = np.append(A1_h1, A1_h2, axis=0)
 
 # Add velocity perturbation
 mask_perturb = np.logical_and(
-                               A2_coords[:, 1] > perturbation_region[0],
-                               A2_coords[:, 1] < perturbation_region[1],
-                              )
-A2_vel[mask_perturb, 1] = dv * (1 + np.cos(8 * np.pi * (A2_coords[mask_perturb, 0] / (boxsize_factor) + 0.25))) * (1 + np.cos(5 * np.pi * (A2_coords[mask_perturb, 1] / (boxsize_factor) - 0.5)))
+    A2_coords[:, 1] > perturbation_region[0],
+    A2_coords[:, 1] < perturbation_region[1],
+)
+A2_vel[mask_perturb, 1] = (
+    dv
+    * (1 + np.cos(8 * np.pi * (A2_coords[mask_perturb, 0] / (boxsize_factor) + 0.25)))
+    * (1 + np.cos(5 * np.pi * (A2_coords[mask_perturb, 1] / (boxsize_factor) - 0.5)))
+)
 
 # Write ICs to file
 with h5py.File(fileOutputName, "w") as f:
     # Header
     grp = f.create_group("/Header")
-    grp.attrs["BoxSize"] = [boxsize_xy[0], boxsize_xy[1] + shift_distance, boxsize_depth]
-    grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+    grp.attrs["BoxSize"] = [
+        boxsize_xy[0],
+        boxsize_xy[1] + shift_distance,
+        boxsize_depth,
+    ]
+    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
@@ -166,29 +204,29 @@ with h5py.File(fileOutputName, "w") as f:
     grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
     grp.attrs["Dimension"] = 3
 
-    #Units
+    # Units
     grp = f.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 = f.create_group("/PartType0")
-    ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+    ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
     ds[()] = A2_coords
-    ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+    ds = grp.create_dataset("Velocities", (numPart, 3), "f")
     ds[()] = A2_vel
-    ds = grp.create_dataset('Masses', (numPart, 1), 'f')
-    ds[()] = A1_m.reshape((numPart,1))
-    ds = grp.create_dataset('Density', (numPart,1), 'f')
-    ds[()] = A1_rho.reshape((numPart,1))
-    ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-    ds[()] = A1_h.reshape((numPart,1))
-    ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-    ds[()] = A1_u.reshape((numPart,1))
-    ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
-    ds[()] = A1_ids.reshape((numPart,1))
-    ds = grp.create_dataset('MaterialIDs', (numPart,1), 'i')
-    ds[()] = A1_mat.reshape((numPart,1))
+    ds = grp.create_dataset("Masses", (numPart, 1), "f")
+    ds[()] = A1_m.reshape((numPart, 1))
+    ds = grp.create_dataset("Density", (numPart, 1), "f")
+    ds[()] = A1_rho.reshape((numPart, 1))
+    ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+    ds[()] = A1_h.reshape((numPart, 1))
+    ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+    ds[()] = A1_u.reshape((numPart, 1))
+    ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+    ds[()] = A1_ids.reshape((numPart, 1))
+    ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
+    ds[()] = A1_mat.reshape((numPart, 1))
diff --git a/examples/Planetary/RayleighTaylor_IdealGas_3D/plotSnapshots.py b/examples/Planetary/RayleighTaylor_IdealGas_3D/plotSnapshots.py
index 99787dea948cdc5f535ccb4d2e1f3d24459ceea9..b1eee9a59b7fefbef582ef8667597d212169ecf5 100644
--- a/examples/Planetary/RayleighTaylor_IdealGas_3D/plotSnapshots.py
+++ b/examples/Planetary/RayleighTaylor_IdealGas_3D/plotSnapshots.py
@@ -27,6 +27,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def make_axes():
     # Use Gridspec to set up figure
     n_gs_ax_x = 40
@@ -50,9 +51,26 @@ def make_axes():
     )
 
     ax_0 = plt.subplot(gs[:n_gs_ax_y, :n_gs_ax_x])
-    ax_1 = plt.subplot(gs[:n_gs_ax_y, n_gs_ax_x+n_gs_ax_gap:2*n_gs_ax_x+n_gs_ax_gap])
-    ax_2 = plt.subplot(gs[:n_gs_ax_y, 2*n_gs_ax_x+2*n_gs_ax_gap:3*n_gs_ax_x+2*n_gs_ax_gap])
-    cax = plt.subplot(gs[:n_gs_ax_y, 3*n_gs_ax_x+2*n_gs_ax_gap+n_gs_cbar_gap:3*n_gs_ax_x+2*n_gs_ax_gap+n_gs_cbar_gap+n_gs_cbar])
+    ax_1 = plt.subplot(
+        gs[:n_gs_ax_y, n_gs_ax_x + n_gs_ax_gap : 2 * n_gs_ax_x + n_gs_ax_gap]
+    )
+    ax_2 = plt.subplot(
+        gs[
+            :n_gs_ax_y,
+            2 * n_gs_ax_x + 2 * n_gs_ax_gap : 3 * n_gs_ax_x + 2 * n_gs_ax_gap,
+        ]
+    )
+    cax = plt.subplot(
+        gs[
+            :n_gs_ax_y,
+            3 * n_gs_ax_x
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax_x
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
 
     axs = [ax_0, ax_1, ax_2]
 
@@ -116,7 +134,6 @@ def plot_kh(ax, snap, cmap, norm):
     ax.set_ylim((0.05 * boxsize_y, 0.95 * boxsize_y))
 
 
-
 if __name__ == "__main__":
 
     # Set colormap
@@ -137,7 +154,14 @@ if __name__ == "__main__":
         time = times[i]
 
         plot_kh(ax, snap, cmap, norm)
-        ax.text(0.5,-0.05,r"$t =\;$"+time, horizontalalignment='center',size=18, transform=ax.transAxes)
+        ax.text(
+            0.5,
+            -0.05,
+            r"$t =\;$" + time,
+            horizontalalignment="center",
+            size=18,
+            transform=ax.transAxes,
+        )
 
     # Colour bar
     sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
diff --git a/examples/Planetary/RayleighTaylor_JupiterLike_3D/makeIC.py b/examples/Planetary/RayleighTaylor_JupiterLike_3D/makeIC.py
index 2aa8a45ea8346bc5d47a6b376beae77e9ec0ee1e..c5e0ee20d130856bf6df1be3c64127734feabb6a 100644
--- a/examples/Planetary/RayleighTaylor_JupiterLike_3D/makeIC.py
+++ b/examples/Planetary/RayleighTaylor_JupiterLike_3D/makeIC.py
@@ -1,59 +1,68 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
- #               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) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+#               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/>.
+#
+##############################################################################
 
 import h5py
 import numpy as np
 import woma
 
 # Load EoS tables
-woma.load_eos_tables(["CD21_HHe","AQUA"])
+woma.load_eos_tables(["CD21_HHe", "AQUA"])
 
 # Generates a swift IC file for the Rayleigh-Taylor instability test
 
 # Constants
-R_earth = 6371000 # Earth radius
-R_jupiter = 11.2089 * R_earth   # Jupiter radius
+R_earth = 6371000  # Earth radius
+R_jupiter = 11.2089 * R_earth  # Jupiter radius
 
 # Parameters
-N2_l = 64             # Number of particles along one edge in lower region
-N2_depth = 18         # Number of particles along in z dimension in lower region
-matID1 = 304          # Upper region material ID: AQUA
-matID2 = 307          # Lower region material ID: CD21 H--He
-rho1_approx = 9000    # Approximate density of upper region. To be recalculated
-rho2 = 3500           # Density of lower region
-g = -31.44            # Constant gravitational acceleration
-P0 = 3.2e+12          # Pressure at interface
+N2_l = 64  # Number of particles along one edge in lower region
+N2_depth = 18  # Number of particles along in z dimension in lower region
+matID1 = 304  # Upper region material ID: AQUA
+matID2 = 307  # Lower region material ID: CD21 H--He
+rho1_approx = 9000  # Approximate density of upper region. To be recalculated
+rho2 = 3500  # Density of lower region
+g = -31.44  # Constant gravitational acceleration
+P0 = 3.2e12  # Pressure at interface
 boxsize_factor = 0.1 * R_jupiter
 dv = 0.00025 * boxsize_factor  # Size of velocity perturbation
-boxsize_xy = [0.5 * boxsize_factor, 1.0 * boxsize_factor]  # Size of the box in x and y dimensions
-boxsize_depth = boxsize_xy[0] * N2_depth / N2_l # Size of simulation box in z dimension
-fixed_region = [0.05 * boxsize_factor, 0.95 * boxsize_factor]       # y-range of non fixed_region particles
-perturbation_region = [0.3 * boxsize_factor, 0.7 * boxsize_factor]  # y-range for the velocity perturbation
+boxsize_xy = [
+    0.5 * boxsize_factor,
+    1.0 * boxsize_factor,
+]  # Size of the box in x and y dimensions
+boxsize_depth = boxsize_xy[0] * N2_depth / N2_l  # Size of simulation box in z dimension
+fixed_region = [
+    0.05 * boxsize_factor,
+    0.95 * boxsize_factor,
+]  # y-range of non fixed_region particles
+perturbation_region = [
+    0.3 * boxsize_factor,
+    0.7 * boxsize_factor,
+]  # y-range for the velocity perturbation
 fileOutputName = "rayleigh_taylor.hdf5"
-#---------------------------------------------------
+# ---------------------------------------------------
 
 # Start by generating grids of particles of the two densities
 numPart2 = N2_l * N2_l * N2_depth
 numPart1 = int(numPart2 / rho2 * rho1_approx)
 N1_l = int(np.cbrt(boxsize_xy[0] * numPart1 / boxsize_depth))
-N1_l -= N1_l % 4 # Make RT symmetric across centre of both instability regions
+N1_l -= N1_l % 4  # Make RT symmetric across centre of both instability regions
 N1_depth = int(boxsize_depth * N1_l / boxsize_xy[0])
 numPart1 = int(N1_l * N1_l * N1_depth)
 numPart = numPart2 + numPart1
@@ -90,13 +99,22 @@ for i in range(N1_depth):
     for j in range(N1_l):
         for k in range(N1_l):
             index = i * N1_l**2 + j * N1_l + k
-            A2_coords1[index, 0] = (j / float(N1_l) + 1. / (2. * N1_l)) * boxsize_xy[0]
-            A2_coords1[index, 1] = (k / float(N1_l) + 1. / (2. * N1_l)) * (0.5 * boxsize_xy[1]) + 0.5 * boxsize_xy[1]
-            A2_coords1[index, 2] = (i / float(N1_depth) + 1. / (2. * N1_depth)) * boxsize_depth
+            A2_coords1[index, 0] = (j / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_xy[
+                0
+            ]
+            A2_coords1[index, 1] = (k / float(N1_l) + 1.0 / (2.0 * N1_l)) * (
+                0.5 * boxsize_xy[1]
+            ) + 0.5 * boxsize_xy[1]
+            A2_coords1[index, 2] = (
+                i / float(N1_depth) + 1.0 / (2.0 * N1_depth)
+            ) * boxsize_depth
             A1_rho1[index] = rho1
 
             # If in top and bottom where particles are fixed
-            if A2_coords1[index, 1] < fixed_region[0] or A2_coords1[index, 1] > fixed_region[1]:
+            if (
+                A2_coords1[index, 1] < fixed_region[0]
+                or A2_coords1[index, 1] > fixed_region[1]
+            ):
                 A1_ids[index] = boundary_particles
                 boundary_particles += 1
 
@@ -106,22 +124,34 @@ for i in range(N2_depth):
     for j in range(N2_l):
         for k in range(N2_l):
             index = i * N2_l**2 + j * N2_l + k
-            A2_coords2[index, 0] = (j / float(N2_l) + 1. / (2. * N2_l)) * boxsize_xy[0]
-            A2_coords2[index, 1] = (k / float(N2_l) + 1. / (2. * N2_l)) * (0.5 * boxsize_xy[1])
-            A2_coords2[index, 2] = (i / float(N2_depth) + 1. / (2. * N2_depth)) * boxsize_depth
+            A2_coords2[index, 0] = (j / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_xy[
+                0
+            ]
+            A2_coords2[index, 1] = (k / float(N2_l) + 1.0 / (2.0 * N2_l)) * (
+                0.5 * boxsize_xy[1]
+            )
+            A2_coords2[index, 2] = (
+                i / float(N2_depth) + 1.0 / (2.0 * N2_depth)
+            ) * boxsize_depth
             A1_rho2[index] = rho2
 
             # If in top and bottom where particles are fixed
-            if A2_coords2[index, 1] < fixed_region[0] or A2_coords2[index, 1] > fixed_region[1]:
+            if (
+                A2_coords2[index, 1] < fixed_region[0]
+                or A2_coords2[index, 1] > fixed_region[1]
+            ):
                 A1_ids[index + numPart1] = boundary_particles
                 boundary_particles += 1
 
 print(
-    "You need to compile the code with " "--enable-boundary-particles=%i" % boundary_particles
+    "You need to compile the code with "
+    "--enable-boundary-particles=%i" % boundary_particles
 )
 
 # Set IDs of non-boundary particles
-A1_ids[A1_ids == 0] = np.linspace(boundary_particles, numPart, numPart - boundary_particles + 1)
+A1_ids[A1_ids == 0] = np.linspace(
+    boundary_particles, numPart, numPart - boundary_particles + 1
+)
 
 # The placement of the lattices are now adjusted to give appropriate interfaces
 # Calculate the separation of particles across the density discontinuity
@@ -130,16 +160,16 @@ pcl_separation_1 = np.cbrt(mass / rho1)
 boundary_separation = 0.5 * (pcl_separation_2 + pcl_separation_1)
 
 # Shift top lattice
-min_y1 = min(A2_coords1[:,1])
+min_y1 = min(A2_coords1[:, 1])
 max_y2 = max(A2_coords2[:, 1])
 shift_distance = boundary_separation - (min_y1 - max_y2)
-A2_coords1[:,1] += shift_distance
+A2_coords1[:, 1] += shift_distance
 
 # Calculate internal energies
 A1_P1 = P0 + g * A1_rho1 * (A2_coords1[:, 1] - 0.5 * boxsize_xy[1])
 A1_P2 = P0 + g * A1_rho2 * (A2_coords2[:, 1] - 0.5 * boxsize_xy[1])
-A1_u1 =  woma.A1_Z_rho_Y(A1_rho1, A1_P1, A1_mat1, Z_choice="u", Y_choice="P")
-A1_u2 =  woma.A1_Z_rho_Y(A1_rho2, A1_P2, A1_mat2, Z_choice="u", Y_choice="P")
+A1_u1 = woma.A1_Z_rho_Y(A1_rho1, A1_P1, A1_mat1, Z_choice="u", Y_choice="P")
+A1_u2 = woma.A1_Z_rho_Y(A1_rho2, A1_P2, A1_mat2, Z_choice="u", Y_choice="P")
 
 # Now the two lattices can be combined
 A2_coords = np.append(A2_coords1, A2_coords2, axis=0)
@@ -152,18 +182,26 @@ A1_h = np.append(A1_h1, A1_h2, axis=0)
 
 # Add velocity perturbation
 mask_perturb = np.logical_and(
-                               A2_coords[:, 1] > perturbation_region[0],
-                               A2_coords[:, 1] < perturbation_region[1],
-                              )
-A2_vel[mask_perturb, 1] = dv * (1 + np.cos(8 * np.pi * (A2_coords[mask_perturb, 0] / (boxsize_factor) + 0.25))) * (1 + np.cos(5 * np.pi * (A2_coords[mask_perturb, 1] / (boxsize_factor) - 0.5)))
+    A2_coords[:, 1] > perturbation_region[0],
+    A2_coords[:, 1] < perturbation_region[1],
+)
+A2_vel[mask_perturb, 1] = (
+    dv
+    * (1 + np.cos(8 * np.pi * (A2_coords[mask_perturb, 0] / (boxsize_factor) + 0.25)))
+    * (1 + np.cos(5 * np.pi * (A2_coords[mask_perturb, 1] / (boxsize_factor) - 0.5)))
+)
 
 
 # Write ICs to file
 with h5py.File(fileOutputName, "w") as f:
     # Header
     grp = f.create_group("/Header")
-    grp.attrs["BoxSize"] = [boxsize_xy[0], boxsize_xy[1] + shift_distance, boxsize_depth]
-    grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+    grp.attrs["BoxSize"] = [
+        boxsize_xy[0],
+        boxsize_xy[1] + shift_distance,
+        boxsize_depth,
+    ]
+    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
@@ -172,29 +210,29 @@ with h5py.File(fileOutputName, "w") as f:
     grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
     grp.attrs["Dimension"] = 3
 
-    #Units
+    # Units
     grp = f.create_group("/Units")
-    grp.attrs["Unit length in cgs (U_L)"] = 100.
-    grp.attrs["Unit mass in cgs (U_M)"] = 1000.
-    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)"] = 100.0
+    grp.attrs["Unit mass in cgs (U_M)"] = 1000.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 = f.create_group("/PartType0")
-    ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+    ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
     ds[()] = A2_coords
-    ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+    ds = grp.create_dataset("Velocities", (numPart, 3), "f")
     ds[()] = A2_vel
-    ds = grp.create_dataset('Masses', (numPart, 1), 'f')
-    ds[()] = A1_m.reshape((numPart,1))
-    ds = grp.create_dataset('Density', (numPart,1), 'f')
-    ds[()] = A1_rho.reshape((numPart,1))
-    ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-    ds[()] = A1_h.reshape((numPart,1))
-    ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-    ds[()] = A1_u.reshape((numPart,1))
-    ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
-    ds[()] = A1_ids.reshape((numPart,1))
-    ds = grp.create_dataset('MaterialIDs', (numPart,1), 'i')
-    ds[()] = A1_mat.reshape((numPart,1))
+    ds = grp.create_dataset("Masses", (numPart, 1), "f")
+    ds[()] = A1_m.reshape((numPart, 1))
+    ds = grp.create_dataset("Density", (numPart, 1), "f")
+    ds[()] = A1_rho.reshape((numPart, 1))
+    ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+    ds[()] = A1_h.reshape((numPart, 1))
+    ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+    ds[()] = A1_u.reshape((numPart, 1))
+    ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+    ds[()] = A1_ids.reshape((numPart, 1))
+    ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
+    ds[()] = A1_mat.reshape((numPart, 1))
diff --git a/examples/Planetary/RayleighTaylor_JupiterLike_3D/plotSnapshots.py b/examples/Planetary/RayleighTaylor_JupiterLike_3D/plotSnapshots.py
index 92ec546b39bb01948bb77c4c3e2e8101556d35c0..da3cd0b147b5b0c52016036e3b6fc6d21af985df 100644
--- a/examples/Planetary/RayleighTaylor_JupiterLike_3D/plotSnapshots.py
+++ b/examples/Planetary/RayleighTaylor_JupiterLike_3D/plotSnapshots.py
@@ -27,6 +27,7 @@ import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
 
+
 def make_axes():
     # Use Gridspec to set up figure
     n_gs_ax_x = 40
@@ -50,11 +51,38 @@ def make_axes():
     )
 
     ax_0 = plt.subplot(gs[:n_gs_ax_y, :n_gs_ax_x])
-    ax_1 = plt.subplot(gs[:n_gs_ax_y, n_gs_ax_x+n_gs_ax_gap:2*n_gs_ax_x+n_gs_ax_gap])
-    ax_2 = plt.subplot(gs[:n_gs_ax_y, 2*n_gs_ax_x+2*n_gs_ax_gap:3*n_gs_ax_x+2*n_gs_ax_gap])
+    ax_1 = plt.subplot(
+        gs[:n_gs_ax_y, n_gs_ax_x + n_gs_ax_gap : 2 * n_gs_ax_x + n_gs_ax_gap]
+    )
+    ax_2 = plt.subplot(
+        gs[
+            :n_gs_ax_y,
+            2 * n_gs_ax_x + 2 * n_gs_ax_gap : 3 * n_gs_ax_x + 2 * n_gs_ax_gap,
+        ]
+    )
 
-    cax_0 = plt.subplot(gs[:int(0.5*n_gs_ax_y)-1, 3*n_gs_ax_x+2*n_gs_ax_gap+n_gs_cbar_gap:3*n_gs_ax_x+2*n_gs_ax_gap+n_gs_cbar_gap+n_gs_cbar])
-    cax_1 = plt.subplot(gs[int(0.5*n_gs_ax_y)+1:, 3*n_gs_ax_x+2*n_gs_ax_gap+n_gs_cbar_gap:3*n_gs_ax_x+2*n_gs_ax_gap+n_gs_cbar_gap+n_gs_cbar])
+    cax_0 = plt.subplot(
+        gs[
+            : int(0.5 * n_gs_ax_y) - 1,
+            3 * n_gs_ax_x
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax_x
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
+    cax_1 = plt.subplot(
+        gs[
+            int(0.5 * n_gs_ax_y) + 1 :,
+            3 * n_gs_ax_x
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax_x
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
 
     axs = [ax_0, ax_1, ax_2]
     caxs = [cax_0, cax_1]
@@ -69,8 +97,8 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
 
     with h5py.File(snap_file, "r") as f:
         # Units from file metadata to SI
-        m=float(f["Units"].attrs["Unit mass in cgs (U_M)"][0]) * 1e-3
-        l=float(f["Units"].attrs["Unit length in cgs (U_L)"][0]) * 1e-2
+        m = float(f["Units"].attrs["Unit mass in cgs (U_M)"][0]) * 1e-3
+        l = float(f["Units"].attrs["Unit length in cgs (U_L)"][0]) * 1e-2
 
         boxsize_x = f["Header"].attrs["BoxSize"][0] * l
         boxsize_y = f["Header"].attrs["BoxSize"][1] * l
@@ -140,21 +168,20 @@ def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
     ax.set_ylim((0.05 * boxsize_y, 0.95 * boxsize_y))
 
 
-
 if __name__ == "__main__":
 
     # Set colormap
-    cmap1 = plt.get_cmap('winter_r')
+    cmap1 = plt.get_cmap("winter_r")
     mat_id1 = 304
     rho_min1 = 6800
     rho_max1 = 8700
-    norm1 = mpl.colors.Normalize(vmin=rho_min1,vmax=rho_max1)
+    norm1 = mpl.colors.Normalize(vmin=rho_min1, vmax=rho_max1)
 
-    cmap2 = plt.get_cmap('autumn')
+    cmap2 = plt.get_cmap("autumn")
     mat_id2 = 307
     rho_min2 = 3450
     rho_max2 = 4050
-    norm2 = mpl.colors.Normalize(vmin=rho_min2,vmax=rho_max2)
+    norm2 = mpl.colors.Normalize(vmin=rho_min2, vmax=rho_max2)
 
     # Generate axes
     axs, caxs = make_axes()
@@ -169,7 +196,14 @@ if __name__ == "__main__":
         time = times[i]
 
         plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2)
-        ax.text(0.5,-0.05,r"$t =\;$"+time+r"$\,$s", horizontalalignment='center',size=18, transform=ax.transAxes)
+        ax.text(
+            0.5,
+            -0.05,
+            r"$t =\;$" + time + r"$\,$s",
+            horizontalalignment="center",
+            size=18,
+            transform=ax.transAxes,
+        )
 
     # Colour bar
     sm1 = plt.cm.ScalarMappable(cmap=cmap1, norm=norm1)
diff --git a/examples/Planetary/SodShock_3D/plotSolution.py b/examples/Planetary/SodShock_3D/plotSolution.py
index 72d16a4b3a4ee2aae6342613317d59085891c9c3..3f95931bd322c3465f272af00bc228fd7d464941 100644
--- a/examples/Planetary/SodShock_3D/plotSolution.py
+++ b/examples/Planetary/SodShock_3D/plotSolution.py
@@ -93,34 +93,34 @@ 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)
+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)
 
 if plot_diffusion:
     alpha_diff_bin, _, _ = stats.binned_statistic(
         x, diffusion, statistic="mean", bins=x_bin_edge
     )
     alpha2_diff_bin, _, _ = stats.binned_statistic(
-        x, diffusion ** 2, statistic="mean", bins=x_bin_edge
+        x, diffusion**2, statistic="mean", bins=x_bin_edge
     )
-    alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin ** 2)
+    alpha_diff_sigma_bin = np.sqrt(alpha2_diff_bin - alpha_diff_bin**2)
 
 if plot_viscosity:
     alpha_visc_bin, _, _ = stats.binned_statistic(
         x, viscosity, statistic="mean", bins=x_bin_edge
     )
     alpha2_visc_bin, _, _ = stats.binned_statistic(
-        x, viscosity ** 2, statistic="mean", bins=x_bin_edge
+        x, viscosity**2, statistic="mean", bins=x_bin_edge
     )
-    alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin ** 2)
+    alpha_visc_sigma_bin = np.sqrt(alpha2_visc_bin - alpha_visc_bin**2)
 
 # Prepare reference solution
 solver = RiemannSolver(gas_gamma)
@@ -131,7 +131,7 @@ rho_s, v_s, P_s, _ = solver.solve(rho_L, v_L, P_L, rho_R, v_R, P_R, x_s / time)
 
 # Additional arrays
 u_s = P_s / (rho_s * (gas_gamma - 1.0))  # internal energy
-s_s = P_s / rho_s ** gas_gamma  # entropic function
+s_s = P_s / rho_s**gas_gamma  # entropic function
 
 # Plot the interesting quantities
 figure(figsize=(7, 7 / 1.6))
diff --git a/examples/Planetary/SquareTest_3D/makeICEqualMass.py b/examples/Planetary/SquareTest_3D/makeICEqualMass.py
index b185edb08d5ce0fb9f502240361f702f326bc7f2..b8e99d621cd36043603769ee7f7361f911765ce8 100644
--- a/examples/Planetary/SquareTest_3D/makeICEqualMass.py
+++ b/examples/Planetary/SquareTest_3D/makeICEqualMass.py
@@ -27,7 +27,7 @@ import numpy as np
 # Parameters
 N_l = 40  # Number of particles on one side
 gamma = 5.0 / 3.0  # Gas adiabatic index
-rho_in_approx = 16**3/10**3  # Density of inner region
+rho_in_approx = 16**3 / 10**3  # Density of inner region
 rho_out = 1.0  # Density of outer region
 P_in = 2.5  # Pressure of inner region
 P_out = 2.5  # Pressure of outer region
diff --git a/examples/Planetary/SquareTest_3D/plotSolution.py b/examples/Planetary/SquareTest_3D/plotSolution.py
index f1010413b758596e9a08769dbdf517655a8c5188..856c2b6b884b8582e30a1b72c8ed3f3cd4ebddbd 100644
--- a/examples/Planetary/SquareTest_3D/plotSolution.py
+++ b/examples/Planetary/SquareTest_3D/plotSolution.py
@@ -150,4 +150,6 @@ if __name__ == "__main__":
     A1_size = size_factor * (A1_m_slice / A1_rho_slice) ** (2 / 3)
 
     # Plot figure
-    plot_square(A1_x_slice, A1_y_slice, A1_rho_slice, A1_u_slice, A1_P_slice, A1_size, boxsize_l)
+    plot_square(
+        A1_x_slice, A1_y_slice, A1_rho_slice, A1_u_slice, A1_P_slice, A1_size, boxsize_l
+    )