examples/HydroTests/PerturbedBox_3D/run.sh

@@ -4,7 +4,7 @@
 if [ ! -e perturbedBox.hdf5 ]
 then
     echo "Generating initial conditions for the perturbed box example..."
-    python makeIC.py 50
+    python3 makeIC.py 50
 fi
 
-../../swift --hydro --threads=16 perturbedBox.yml 2>&1 | tee output.log
+../../../swift --hydro --threads=16 perturbedBox.yml 2>&1 | tee output.log

examples/HydroTests/Rayleigh-Taylor_2D/makeIC.py

@@ -26,12 +26,12 @@ from scipy.optimize import bisect
 
 # Parameters
 N = [128, 192]  # Particles along one edge
-gamma = 7./5.  # Gas adiabatic index
-dv = 0.025    # velocity perturbation
-rho_h = 2    # high density region
-rho_l = 1    # Low density region
-g = -0.5    # gravitational acceleration
-box_size = [1., 1.5]  # size of the box
+gamma = 7.0 / 5.0  # Gas adiabatic index
+dv = 0.025  # velocity perturbation
+rho_h = 2  # high density region
+rho_l = 1  # Low density region
+g = -0.5  # gravitational acceleration
+box_size = [1.0, 1.5]  # size of the box
 
 fixed = [0.1, 1.4]  # y-range of non fixed particles
 perturbation_box = [0.3, 1.2]  # y-range for the velocity perturbation
@@ -40,7 +40,7 @@ fileOutputName = "rayleigh_taylor.hdf5"
 
 # ---------------------------------------------------
 
-if (N[0] / box_size[0] != N[1] / box_size[1]):
+if N[0] / box_size[0] != N[1] / box_size[1]:
     raise Exception("Suppose the same ratio for each directions")
 
 
@@ -53,11 +53,11 @@ def density(y):
 
     ind = y < 0.5 * box_size[1]
     rho = np.zeros(y.shape)
-    tmp = (gamma - 1.) * g / (gamma * P0)
-    alpha = 1. / (gamma - 1.)
+    tmp = (gamma - 1.0) * g / (gamma * P0)
+    alpha = 1.0 / (gamma - 1.0)
 
-    rho[ind] = rho_l * (1 + rho_l * tmp * (y[ind] - 0.5 * box_size[1]))**alpha
-    rho[~ind] = rho_h * (1 + rho_h * tmp * (y[~ind] - 0.5 * box_size[1]))**alpha
+    rho[ind] = rho_l * (1 + rho_l * tmp * (y[ind] - 0.5 * box_size[1])) ** alpha
+    rho[~ind] = rho_h * (1 + rho_h * tmp * (y[~ind] - 0.5 * box_size[1])) ** alpha
 
     return rho
 
@@ -72,14 +72,14 @@ def mass(y):
     ind = y < 0.5 * box_size[1]
     m = np.zeros(y.shape)
 
-    B = (gamma - 1.) * g / (gamma * P0)
-    alpha = 1. / (gamma - 1.)
+    B = (gamma - 1.0) * g / (gamma * P0)
+    alpha = 1.0 / (gamma - 1.0)
 
-    m[ind] = (1 + B * rho_l * (y[ind] - 0.5 * box_size[1]))**(alpha + 1)
+    m[ind] = (1 + B * rho_l * (y[ind] - 0.5 * box_size[1])) ** (alpha + 1)
 
-    m[~ind] = (1 + B * rho_h * (y[~ind] - 0.5 * box_size[1]))**(alpha + 1)
+    m[~ind] = (1 + B * rho_h * (y[~ind] - 0.5 * box_size[1])) ** (alpha + 1)
 
-    m -= (1 - 0.5 * B * box_size[1] * rho_l)**(alpha + 1)
+    m -= (1 - 0.5 * B * box_size[1] * rho_l) ** (alpha + 1)
 
     return box_size[0] * m / (B * (alpha + 1))
 
@@ -115,8 +115,8 @@ def entropy(y):
     ind = y < 0.5 * box_size[1]
 
     a = np.zeros(y.shape)
-    a[ind] = P0 * rho_l**(-gamma)
-    a[~ind] = P0 * rho_h**(-gamma)
+    a[ind] = P0 * rho_l ** (-gamma)
+    a[~ind] = P0 * rho_h ** (-gamma)
     return a
 
 
@@ -132,11 +132,11 @@ def growth_rate():
     Compute the growth rate of the instability.
     Assumes a wavelength equal to the boxsize.
     """
-    ymin = 0.
+    ymin = 0.0
     ymax = box_size[1]
     A = density(ymax) - density(ymin)
     A /= density(ymax) + density(ymin)
-    return np.sqrt(A * np.abs(g) * ymax / (2. * np.pi))
+    return np.sqrt(A * np.abs(g) * ymax / (2.0 * np.pi))
 
 
 def vy(x, y):
@@ -186,20 +186,21 @@ if __name__ == "__main__":
 
             coords[index, 0] = x
             coords[index, 1] = y_j
-            if (y_j < fixed[0] or y_j > fixed[1]):
+            if y_j < fixed[0] or y_j > fixed[1]:
                 ids[index] = uni_id
                 uni_id += 1
 
-    print("You need to compile the code with "
-          "--enable-boundary-particles=%i" % uni_id)
-    ids[ids == 0] = np.linspace(uni_id, numPart, numPart-uni_id+1)
+    print(
+        "You need to compile the code with " "--enable-boundary-particles=%i" % uni_id
+    )
+    ids[ids == 0] = np.linspace(uni_id, numPart, numPart - uni_id + 1)
 
     # density
     rho = density(coords[:, 1])
 
     # internal energy
     a = entropy(coords[:, 1])
-    u = a * rho**(gamma-1) / (gamma - 1.)
+    u = a * rho ** (gamma - 1) / (gamma - 1.0)
 
     # smoothing length
     h = smoothing_length(rho, m)
@@ -208,7 +209,7 @@ if __name__ == "__main__":
     vel[:, 1] = vy(coords[:, 0], coords[:, 1])
 
     # File
-    fileOutput = h5py.File(fileOutputName, 'w')
+    fileOutput = h5py.File(fileOutputName, "w")
 
     # Header
     grp = fileOutput.create_group("/Header")
@@ -224,27 +225,27 @@ if __name__ == "__main__":
 
     # Units
     grp = fileOutput.create_group("/Units")
-    grp.attrs["Unit length in cgs (U_L)"] = 1.
-    grp.attrs["Unit mass in cgs (U_M)"] = 1.
-    grp.attrs["Unit time in cgs (U_t)"] = 1.
-    grp.attrs["Unit current in cgs (U_I)"] = 1.
-    grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+    grp.attrs["Unit length in cgs (U_L)"] = 1.0
+    grp.attrs["Unit mass in cgs (U_M)"] = 1.0
+    grp.attrs["Unit time in cgs (U_t)"] = 1.0
+    grp.attrs["Unit current in cgs (U_I)"] = 1.0
+    grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
 
     # Particle group
     grp = fileOutput.create_group("/PartType0")
-    ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+    ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
     ds[()] = coords
-    ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+    ds = grp.create_dataset("Velocities", (numPart, 3), "f")
     ds[()] = vel
-    ds = grp.create_dataset('Masses', (numPart, 1), 'f')
+    ds = grp.create_dataset("Masses", (numPart, 1), "f")
     ds[()] = m.reshape((numPart, 1))
-    ds = grp.create_dataset('SmoothingLength', (numPart, 1), 'f')
+    ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
     ds[()] = h.reshape((numPart, 1))
-    ds = grp.create_dataset('InternalEnergy', (numPart, 1), 'f')
+    ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
     ds[()] = u.reshape((numPart, 1))
-    ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
+    ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
     ds[()] = ids.reshape((numPart, 1))
-    ds = grp.create_dataset('Density', (numPart, 1), 'f')
+    ds = grp.create_dataset("Density", (numPart, 1), "f")
     ds[()] = rho.reshape((numPart, 1))
 
     fileOutput.close()

examples/HydroTests/Rayleigh-Taylor_2D/makeMovie.py

@@ -45,6 +45,7 @@ class Metadata(object):
     """
     Copy the useful data in order to decrease the memory usage
     """
+
     def __init__(self, data):
         metadata = data.metadata
         self.t = metadata.t
@@ -91,7 +92,7 @@ def load_and_make_image(filename, res, property):
     border = int(0.2 * m_res)
 
     # first part of the image
-    ylim = np.array([0., 1.])
+    ylim = np.array([0.0, 1.0])
 
     data = load(filename)
     image[:m_res, :m_res] = project(data, m_res, property, ylim)
@@ -144,27 +145,36 @@ def create_movie(filename, start, stop, resolution, property, output_filename):
 
     norm = LogNorm(vmin=vmin, vmax=vmax, clip="black")
 
-    image = ax.imshow(np.zeros_like(frames[0]), origin="lower",
-                      norm=norm)
+    image = ax.imshow(np.zeros_like(frames[0]), origin="lower", norm=norm)
 
     description_text = ax.text(
-        0.5, 0.5,
+        0.5,
+        0.5,
         get_simulation_information(metadata[0]),
-        va="center", ha="center",
-        **text_args, transform=ax.transAxes,
+        va="center",
+        ha="center",
+        **text_args,
+        transform=ax.transAxes,
     )
 
     time_text = ax.text(
-        (1 - 0.025 * 0.25), 0.975,
+        (1 - 0.025 * 0.25),
+        0.975,
         time_formatter(metadata[0]),
         **text_args,
-        va="top", ha="right",
+        va="top",
+        ha="right",
         transform=ax.transAxes,
     )
 
     ax.text(
-        0.025 * 0.25, 0.975, name, **text_args, va="top", ha="left",
-        transform=ax.transAxes
+        0.025 * 0.25,
+        0.975,
+        name,
+        **text_args,
+        va="top",
+        ha="left",
+        transform=ax.transAxes,
     )
 
     def frame(n):
@@ -174,19 +184,19 @@ def create_movie(filename, start, stop, resolution, property, output_filename):
             time_text.set_text(time_formatter(metadata[n]))
 
         if generate_png:
-            name = filename + "_{:04d}".format(n+info_frames)
+            name = filename + "_{:04d}".format(n + info_frames)
             fig.savefig(name + ".png")
 
         else:
             return (image,)
 
     if generate_png:
-        for i in range(-info_frames, stop-start):
+        for i in range(-info_frames, stop - start):
             frame(i)
     else:
-        animation = FuncAnimation(fig, frame,
-                                  range(-info_frames, stop-start),
-                                  interval=40)
+        animation = FuncAnimation(
+            fig, frame, range(-info_frames, stop - start), interval=40
+        )
         animation.save(output_filename)
 
 
@@ -232,13 +242,11 @@ if __name__ == "__main__":
     )
 
     parser.add_argument(
-        "-i", "--initial", help="Initial snapshot. Default: 0",
-        type=int, default=0
+        "-i", "--initial", help="Initial snapshot. Default: 0", type=int, default=0
     )
 
     parser.add_argument(
-        "-f", "--final", help="Final snapshot. Default: 40",
-        type=int, default=40
+        "-f", "--final", help="Final snapshot. Default: 40", type=int, default=40
     )
 
     parser.add_argument(

examples/HydroTests/Rayleigh-Taylor_2D/run.sh

@@ -6,10 +6,10 @@ set -e
 if [ ! -e rayleigh_taylor.hdf5 ]
 then
     echo "Generating initial conditions for the Rayleigh Taylor example..."
-    python makeIC.py
+    python3 makeIC.py
 fi
 
 # Run SWIFT
-../../swift --hydro --external-gravity  --threads=8 rayleigh_taylor.yml 2>&1 | tee output.log
+../../../swift --hydro --external-gravity  --threads=8 rayleigh_taylor.yml 2>&1 | tee output.log
 
-python makeMovie.py -i 0 -f 1001
+python3 makeMovie.py -i 0 -f 1001