examples/HydroTests/KeplerianRing/plotSolution.py

@@ -26,11 +26,13 @@
 # and then our own density as a function of radius calculation.
 
 import matplotlib
+
 matplotlib.use("Agg")
 matplotlib.rc("text", usetex=True)
 
 try:
     import yt
+
     ytavail = True
 except ImportError:
     print("yt not found. Falling back on homebrew plots.")
@@ -66,7 +68,7 @@ def get_axes_grid(figure):
     gs = gridspec.GridSpec(2, 30)
 
     grid = []
-    
+
     grid.append(figure.add_subplot(gs[0, 0:10]))
     grid.append(figure.add_subplot(gs[0, 10:20]))
     grid.append(figure.add_subplot(gs[0, 20:]))
@@ -103,7 +105,7 @@ def chi_square(observed, expected):
     masked_expected = np.array(expected)[mask]
     masked_observed = np.array(observed)[mask]
 
-    return sum(((masked_observed - masked_expected)**2)/masked_expected**2)
+    return sum(((masked_observed - masked_expected) ** 2) / masked_expected ** 2)
 
 
 def load_data(filename):
@@ -118,11 +120,13 @@ def load_data(filename):
         # Check if z = 0 for all particles. If so we need to set the cetnre
         # to have z = 0 also.
         if np.sum(file["PartType0"]["Coordinates"][:, 2]) == 0:
-            centre = [boxsize[0] / 2., boxsize[0] / 2., 0.]
+            centre = [boxsize[0] / 2.0, boxsize[0] / 2.0, 0.0]
         else:
-            centre = boxsize / 2.
+            centre = boxsize / 2.0
 
-        radii = np.sqrt(np.sum(((file["PartType0"]["Coordinates"][...] - centre).T)**2, 0))
+        radii = np.sqrt(
+            np.sum(((file["PartType0"]["Coordinates"][...] - centre).T) ** 2, 0)
+        )
         masses = file["PartType0"]["Masses"][...]
 
     return radii, masses
@@ -137,16 +141,16 @@ def bin_density_r(radii, density, binrange, binnumber):
     indicies = np.digitize(radii, bins)
 
     binned_masses = np.zeros(len(bins) - 1)
-    
+
     for index, bin in enumerate(indicies):
         if bin >= len(bins) - 1:
             continue
 
         binned_masses[bin] += density[index]
 
-    areas = [np.pi * (a**2 - b**2) for a, b in zip(bins[1:], bins[:-1])]
-    binned_densities = binned_masses/areas
-    
+    areas = [np.pi * (a ** 2 - b ** 2) for a, b in zip(bins[1:], bins[:-1])]
+    binned_densities = binned_masses / areas
+
     return bins, binned_densities
 
 
@@ -179,7 +183,9 @@ def get_mass_outside_inside(snap, radius=2, filename="keplerian_ring"):
     return below, above
 
 
-def get_derived_data(minsnap, maxsnap, filename="keplerian_ring", binnumber=50, radius=2):
+def get_derived_data(
+    minsnap, maxsnap, filename="keplerian_ring", binnumber=50, radius=2
+):
     """
     Gets the derived data from our snapshots, i.e. the
     density(r) profile and the chi squared (based on the
@@ -188,36 +194,32 @@ def get_derived_data(minsnap, maxsnap, filename="keplerian_ring", binnumber=50,
 
     initial = get_density_r(minsnap, filename, binnumber=binnumber)
     other_densities = [
-        get_density_r(snap, binnumber=binnumber)[1] for snap in tqdm(
-            range(minsnap+1, maxsnap+1), desc="Densities"
-        )
+        get_density_r(snap, binnumber=binnumber)[1]
+        for snap in tqdm(range(minsnap + 1, maxsnap + 1), desc="Densities")
     ]
     densities = [initial[1]] + other_densities
 
     masses_inside_outside = [
-        get_mass_outside_inside(snap, radius=radius, filename=filename)[0] for snap in tqdm(
-            range(minsnap, maxsnap+1), desc="Mass Flow"
-        )
+        get_mass_outside_inside(snap, radius=radius, filename=filename)[0]
+        for snap in tqdm(range(minsnap, maxsnap + 1), desc="Mass Flow")
     ]
 
     # Between the initial conditions and the first snapshot we hope that there
-    # has been no mass flow, hence the [0.] + 
-    mass_flows = [0.] + [
-        y - x for x, y in zip(
-            masses_inside_outside[1:],
-            masses_inside_outside[:-1]
-        )
+    # has been no mass flow, hence the [0.] +
+    mass_flows = [0.0] + [
+        y - x for x, y in zip(masses_inside_outside[1:], masses_inside_outside[:-1])
     ]
 
     cumulative_mass_flows = [
-        sum(mass_flows[:x])/masses_inside_outside[0] for x in range(len(mass_flows))
+        sum(mass_flows[:x]) / masses_inside_outside[0] for x in range(len(mass_flows))
     ]
 
     chisq = [
         chi_square(
-            dens[int(0.1*binnumber):int(0.4*binnumber)],
-            initial[1][int(0.1*binnumber):int(0.4*binnumber)]
-        ) for dens in tqdm(densities, desc="Chi Squared")
+            dens[int(0.1 * binnumber) : int(0.4 * binnumber)],
+            initial[1][int(0.1 * binnumber) : int(0.4 * binnumber)],
+        )
+        for dens in tqdm(densities, desc="Chi Squared")
     ]
 
     return initial[0], densities, chisq, cumulative_mass_flows
@@ -228,8 +230,10 @@ def plot_chisq(ax, minsnap, maxsnap, chisq, filename="keplerian_ring"):
     Plot the chisq(rotation).
     """
     snapshots = np.arange(minsnap, maxsnap + 1)
-    rotations = [convert_snapshot_number_to_rotations_at(1, snap, filename) for snap in snapshots]
-    ax.plot(rotations, np.array(chisq)/max(chisq))
+    rotations = [
+        convert_snapshot_number_to_rotations_at(1, snap, filename) for snap in snapshots
+    ]
+    ax.plot(rotations, np.array(chisq) / max(chisq))
 
     ax.set_xlabel("Number of rotations")
     ax.set_ylabel("$\chi^2 / \chi^2_{{max}}$ = {:3.5f}".format(max(chisq)))
@@ -242,7 +246,9 @@ def plot_mass_flow(ax, minsnap, maxsnap, mass_flow, filename="keplerian_ring"):
     Plot the mass_flow(rotation).
     """
     snapshots = np.arange(minsnap, maxsnap + 1)
-    rotations = [convert_snapshot_number_to_rotations_at(1, snap, filename) for snap in snapshots]
+    rotations = [
+        convert_snapshot_number_to_rotations_at(1, snap, filename) for snap in snapshots
+    ]
 
     ax.plot(rotations, mass_flow)
 
@@ -259,7 +265,7 @@ def plot_density_r(ax, bins, densities, snaplist, filename="keplerian_ring"):
     Densities is the _full_ list of density profiles, and
     snaplist is the ones that you wish to plot.
     """
-    radii = [(x + y)/2 for x, y in zip(bins[1:], bins[:-1])]
+    radii = [(x + y) / 2 for x, y in zip(bins[1:], bins[:-1])]
 
     for snap in snaplist:
         index = snap - snaplist[0]
@@ -281,25 +287,34 @@ def plot_extra_info(ax, filename):
     """
 
     metadata = get_metadata(filename)
-    
-    git = metadata['code']['Git Revision'].decode("utf-8")
-    compiler_name = metadata['code']['Compiler Name'].decode("utf-8")
-    compiler_version = metadata['code']['Compiler Version'].decode("utf-8")
-    scheme = metadata['hydro']['Scheme'].decode("utf-8")
-    kernel = metadata['hydro']['Kernel function'].decode("utf-8")
+
+    git = metadata["code"]["Git Revision"].decode("utf-8")
+    compiler_name = metadata["code"]["Compiler Name"].decode("utf-8")
+    compiler_version = metadata["code"]["Compiler Version"].decode("utf-8")
+    scheme = metadata["hydro"]["Scheme"].decode("utf-8")
+    kernel = metadata["hydro"]["Kernel function"].decode("utf-8")
     gas_gamma = metadata["hydro"]["Adiabatic index"][0]
     neighbors = metadata["hydro"]["Kernel target N_ngb"][0]
     eta = metadata["hydro"]["Kernel eta"][0]
-    
 
-    ax.text(-0.49, 0.9, "Keplerian Ring with  $\\gamma={:4.4f}$ in 2/3D".format(gas_gamma), fontsize=11)
+    ax.text(
+        -0.49,
+        0.9,
+        "Keplerian Ring with  $\\gamma={:4.4f}$ in 2/3D".format(gas_gamma),
+        fontsize=11,
+    )
     ax.text(-0.49, 0.8, f"Compiler: {compiler_name} {compiler_version}", fontsize=10)
     ax.text(-0.49, 0.7, "Rotations are quoted at $r=1$", fontsize=10)
-    ax.plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
-    ax.text(-0.49, 0.5, f"$\\textsc{{Swift}}$ {git}", fontsize=10)
+    ax.plot([-0.49, 0.1], [0.62, 0.62], "k-", lw=1)
+    ax.text(-0.49, 0.5, f"$SWIFT$ {git}", fontsize=10)
     ax.text(-0.49, 0.4, scheme, fontsize=10)
     ax.text(-0.49, 0.3, kernel, fontsize=10)
-    ax.text(-0.49, 0.2, "${:2.2f}$ neighbours ($\\eta={:3.3f}$)".format(neighbors, eta), fontsize=10)
+    ax.text(
+        -0.49,
+        0.2,
+        "${:2.2f}$ neighbours ($\\eta={:3.3f}$)".format(neighbors, eta),
+        fontsize=10,
+    )
 
     ax.set_axis_off()
     ax.set_xlim(-0.5, 0.5)
@@ -308,7 +323,9 @@ def plot_extra_info(ax, filename):
     return
 
 
-def surface_density_plot_no_yt(ax, snapnum, filename="keplerian_ring", density_limits=None, vlim=None):
+def surface_density_plot_no_yt(
+    ax, snapnum, filename="keplerian_ring", density_limits=None, vlim=None
+):
     """
     Make the surface density plot (sans yt).
 
@@ -328,7 +345,6 @@ def surface_density_plot_no_yt(ax, snapnum, filename="keplerian_ring", density_l
 
     ax.scatter(x, y, c=density, vmin=vlim[0], vmax=vlim[1], s=0.1)
 
-
     metadata = get_metadata("{}_{:04d}.hdf5".format(filename, snapnum))
     period = metadata["period"]
 
@@ -336,10 +352,9 @@ def surface_density_plot_no_yt(ax, snapnum, filename="keplerian_ring", density_l
         2.5,
         7.5,
         "Snapshot = {:04d}\nRotations = {:1.2f}".format(
-            snapnum,
-            float(metadata["header"]["Time"])/period
+            snapnum, float(metadata["header"]["Time"]) / period
         ),
-        color='black'
+        color="black",
     )
 
     t.set_bbox(dict(alpha=0.5, color="white"))
@@ -351,23 +366,24 @@ def surface_density_plot_no_yt(ax, snapnum, filename="keplerian_ring", density_l
 
     # We now want to remove all of the ticklabels.
 
-    for axis in ['x', 'y']:
+    for axis in ["x", "y"]:
         ax.tick_params(
-            axis=axis,          
-            which='both',      
-            bottom='off',      
-            top='off',         
-            left='off',
-            right='off',
-            labelleft='off',
-            labelbottom='off'
-        ) 
+            axis=axis,
+            which="both",
+            bottom="off",
+            top="off",
+            left="off",
+            right="off",
+            labelleft="off",
+            labelbottom="off",
+        )
 
     return density_limits, vlim
 
 
-
-def surface_density_plot(ax, snapnum, filename="keplerian_ring", density_limits=None, vlim=None):
+def surface_density_plot(
+    ax, snapnum, filename="keplerian_ring", density_limits=None, vlim=None
+):
     """
     Make the surface density plot (via yt).
 
@@ -376,11 +392,7 @@ def surface_density_plot(ax, snapnum, filename="keplerian_ring", density_limits=
     max/min.
     """
 
-    unit_base = {
-        'length': (1.0, 'cm'),
-        'velocity': (1.0, 'cm/s'),
-        'mass': (1.0, 'g')
-    }
+    unit_base = {"length": (1.0, "cm"), "velocity": (1.0, "cm/s"), "mass": (1.0, "g")}
 
     filename = "{}_{:04d}.hdf5".format(filename, snapnum)
 
@@ -391,16 +403,11 @@ def surface_density_plot(ax, snapnum, filename="keplerian_ring", density_limits=
         # number. Just return what we're given.
         return density_limits, vlim
 
-    projection_plot = yt.ProjectionPlot(
-        snap,
-        "z",
-        ("gas", "cell_mass"),
-        width=5.5
-    )
+    projection_plot = yt.ProjectionPlot(snap, "z", ("gas", "cell_mass"), width=5.5)
 
     max_density = snap.all_data()[("gas", "cell_mass")].max()
     min_density = snap.all_data()[("gas", "cell_mass")].min()
-    
+
     new_density_limits = (min_density, max_density)
 
     if density_limits is None:
@@ -417,12 +424,7 @@ def surface_density_plot(ax, snapnum, filename="keplerian_ring", density_limits=
     if vlim is None:
         vlim = new_vlim
 
-    ax.imshow(
-        data[0],
-        cmap=data[1],
-        vmin=vlim[0],
-        vmax=vlim[1]
-    )
+    ax.imshow(data[0], cmap=data[1], vmin=vlim[0], vmax=vlim[1])
 
     metadata = get_metadata(filename)
     period = metadata["period"]
@@ -431,25 +433,24 @@ def surface_density_plot(ax, snapnum, filename="keplerian_ring", density_limits=
         20,
         80,
         "Snapshot = {:04d}\nRotations = {:1.2f}".format(
-            snapnum,
-            float(snap.current_time)/period
+            snapnum, float(snap.current_time) / period
         ),
-        color='white'
+        color="white",
     )
 
     # We now want to remove all of the ticklabels.
 
-    for axis in ['x', 'y']:
+    for axis in ["x", "y"]:
         ax.tick_params(
-            axis=axis,          
-            which='both',      
-            bottom='off',      
-            top='off',         
-            left='off',
-            right='off',
-            labelleft='off',
-            labelbottom='off'
-        ) 
+            axis=axis,
+            which="both",
+            bottom="off",
+            top="off",
+            left="off",
+            right="off",
+            labelleft="off",
+            labelbottom="off",
+        )
 
     return density_limits, vlim
 
@@ -485,7 +486,7 @@ if __name__ == "__main__":
              Default: First snapshot in the folder.
              """,
         default=-1,
-        required=False
+        required=False,
     )
 
     parser.add_argument(
@@ -496,7 +497,7 @@ if __name__ == "__main__":
              Default: (end - beginning) // 2.
              """,
         default=-1,
-        required=False
+        required=False,
     )
 
     parser.add_argument(
@@ -507,7 +508,7 @@ if __name__ == "__main__":
              Default: Last snapshot in the folder.
              """,
         default=-1,
-        required=False
+        required=False,
     )
 
     parser.add_argument(
@@ -518,7 +519,7 @@ if __name__ == "__main__":
              Default: plot.png
              """,
         default="plot.png",
-        required=False
+        required=False,
     )
 
     parser.add_argument(
@@ -529,7 +530,7 @@ if __name__ == "__main__":
              Default: 100.
              """,
         default=100,
-        required=False
+        required=False,
     )
 
     parser.add_argument(
@@ -541,7 +542,7 @@ if __name__ == "__main__":
              Default: 0.
              """,
         default=0,
-        required=False
+        required=False,
     )
 
     parser.add_argument(
@@ -553,9 +554,9 @@ if __name__ == "__main__":
              default value of keplerian_ring.
              """,
         default="keplerian_ring",
-        required=False
+        required=False,
     )
-    
+
     parser.add_argument(
         "-y",
         "--yt",
@@ -565,7 +566,7 @@ if __name__ == "__main__":
              setup. Default: False
              """,
         default=False,
-        required=False
+        required=False,
     )
 
     args = vars(parser.parse_args())
@@ -575,9 +576,9 @@ if __name__ == "__main__":
     if args["beginning"] == -1:
         # We look for the maximum number of snapshots.
         numbers = [
-            int(x[len(filename)+1:-5])
-            for x in os.listdir() if
-            x[:len(filename)] == filename and x[-1] == "5"
+            int(x[len(filename) + 1 : -5])
+            for x in os.listdir()
+            if x[: len(filename)] == filename and x[-1] == "5"
         ]
 
         snapshots = [min(numbers), (max(numbers) - min(numbers)) // 2, max(numbers)]
@@ -593,25 +594,21 @@ if __name__ == "__main__":
     for snap, ax in zip(snapshots, tqdm(axes[0:3], desc="Images")):
         if args["yt"] and ytavail:
             density_limits, vlim = surface_density_plot(
-                ax,
-                snap,
-                density_limits=density_limits,
-                vlim=vlim
+                ax, snap, density_limits=density_limits, vlim=vlim
             )
 
             figure.subplots_adjust(hspace=0, wspace=0)
         else:
             density_limits, vlim = surface_density_plot_no_yt(
-                ax,
-                snap,
-                density_limits=density_limits,
-                vlim=vlim
+                ax, snap, density_limits=density_limits, vlim=vlim
             )
 
     # Now we need to do the density(r) plot.
 
     # Derived data includes density profiles and chi squared
-    derived_data = get_derived_data(snapshots[0], snapshots[2], binnumber=int(args["nbins"]))
+    derived_data = get_derived_data(
+        snapshots[0], snapshots[2], binnumber=int(args["nbins"])
+    )
 
     if args["plotmassflow"]:
         plot_mass_flow(axes[3], snapshots[0], snapshots[2], derived_data[3])
@@ -623,5 +620,3 @@ if __name__ == "__main__":
     plot_extra_info(axes[5], "keplerian_ring_0000.hdf5")
 
     figure.savefig(args["filename"], dpi=300)
-
-

examples/HydroTests/KeplerianRing/run.sh

@@ -9,7 +9,7 @@ then
 fi
 
 rm -rf keplerian_ring_*.hdf5
-../../swift --external-gravity --hydro --threads=1 --verbose=1 keplerian_ring.yml 2>&1 | tee output.log
+../../../swift --external-gravity --hydro --threads=1 --verbose=1 keplerian_ring.yml 2>&1 | tee output.log
 
 echo
 echo

examples/HydroTests/KeplerianRing/testplots.py

@@ -34,12 +34,7 @@ if __name__ == "__main__":
     data = yt.load("initial_conditions.hdf5")
 
     plot = yt.ParticlePlot(
-        data,
-        "particle_position_x",
-        "particle_position_y",
-        "density"
+        data, "particle_position_x", "particle_position_y", "density"
     )
 
     plot.save("test.png")
-
-

examples/HydroTests/KeplerianRing/write_gadget.py

@@ -116,12 +116,12 @@ def write_header(f, boxsize, flag_entropy, np_total, np_total_hw, other=False):
     # We'll first build a dictionary to iterate through.
 
     default_attributes = {
-        "BoxSize" : boxsize,
-        "Flag_Entropy_ICs" : flag_entropy,
-        "NumPart_Total" : np_total,
-        "NumPart_Total_HighWord" : np_total_hw,
-        "NumFilesPerSnapshot" : 1,  # required for backwards compatibility
-        "NumPart_ThisFile" : np_total, # Also required for bw compatibility
+        "BoxSize": boxsize,
+        "Flag_Entropy_ICs": flag_entropy,
+        "NumPart_Total": np_total,
+        "NumPart_Total_HighWord": np_total_hw,
+        "NumFilesPerSnapshot": 1,  # required for backwards compatibility
+        "NumPart_ThisFile": np_total,  # Also required for bw compatibility
     }
 
     if other:
@@ -164,9 +164,7 @@ def write_runtime_pars(f, periodic_boundary, other=False):
 
     # First build the dictionary
 
-    default_attributes = {
-        "PeriodicBoundariesOn" : periodic_boundary,
-    }
+    default_attributes = {"PeriodicBoundariesOn": periodic_boundary}
 
     if other:
         attributes = dict(default_attributes, **other)
@@ -283,16 +281,16 @@ def write_block(f, part_type, pos, vel, ids, mass, int_energy, smoothing, other=
               the particle data. They will be passed such that the key is
               the name of the dataset in the hdf5 file.
     """
-    
+
     # Build the dictionary
 
     default_data = {
-        "Coordinates" : pos,
-        "Velocities" : vel,
-        "ParticleIDs" : ids,
-        "Masses" : mass,
-        "InternalEnergy" : int_energy,
-        "SmoothingLength" : smoothing,
+        "Coordinates": pos,
+        "Velocities": vel,
+        "ParticleIDs": ids,
+        "Masses": mass,
+        "InternalEnergy": int_energy,
+        "SmoothingLength": smoothing,
     }
 
     if other:
@@ -306,4 +304,3 @@ def write_block(f, part_type, pos, vel, ids, mass, int_energy, smoothing, other=
         particles.create_dataset(name, data=value)
 
     return
-

examples/HydroTests/NFW_Hydrostatic/NFW_Hydrostatic.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e21 # kpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
+
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.0  # The end time of the simulation (in internal units).
+  dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
+  dt_max:     0.1  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            snapshot # Common part of the name of output files
+  time_first:          0.       # Time of the first output (in internal units)
+  delta_time:          2e-2     # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-3 # Time between statistics output
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                       0.05    # Constant dimensionless multiplier for time integration.
+  MAC:                       geometric
+  theta_cr:                  0.7     
+  comoving_DM_softening:     0.01 # Comoving softening length (in internal units).
+  max_physical_DM_softening: 0.01    # Physical softening length (in internal units).
+  comoving_baryon_softening:     0.01 # Comoving softening length (in internal units).
+  max_physical_baryon_softening: 0.01    # Physical softening length (in internal units).
+  
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./nfw.hdf5     # The file to read
+  periodic:   0                     # Non-periodic BCs
+  shift:    [0,0,0]   # Centre the box
+
+NFWPotential:
+  useabspos:       0              # 0 -> positions based on centre, 1 -> absolute positions
+  position:        [0.,0.,0.]     # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
+  M_200:            9.5          # M200 of the galaxy disk
+  h:               0.72          # reduced Hubble constant (value does not specify the used units!)
+  concentration:   17.0            # concentration of the Halo
+  diskfraction:    0.15          # Disk mass fraction (here this is the gas)
+  bulgefraction:   0.0            # Bulge mass fraction
+  timestep_mult:   0.01           # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
+  epsilon:         0.001            # Softening size (internal units)

examples/HydroTests/NFW_Hydrostatic/README

+
+Evolve gas at the static equilibrium in an NFW potential for several dynamical times.
+This test allows to compare the impact of different hydro solvers on the initial 
+hydrostatic equilibrium.
+
+To run SWIFT configured with the sphenix hydro solver:
+  ./configure --with-ext-potential=nfw --with-hydro=sphenix  
+  
+To run SWIFT configured with the gadget2 hydro solver:
+  ./configure --with-ext-potential=nfw --with-hydro=gadget2 --disable-hand-vec   
+
+To run SWIFT configured with the gizmo (MFV) hydro solver:
+  ./configure --with-ext-potential=nfw --with-hydro=gizmo-mfv --with-riemann-solver=exact    

examples/HydroTests/NFW_Hydrostatic/getICs.sh

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

examples/HydroTests/NFW_Hydrostatic/plotGasDensityProfile.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2023 Yves Revaz (yves.revaz@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+from glob import glob
+import h5py
+
+
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
+
+MyrInSec = 31557600000000.0
+gcm3InAcc = 1 / 5.978637406556783e23
+
+
+def ComputeDensity(snap):
+
+    # Read the initial state of the gas
+    f = h5py.File(snap, "r")
+
+    # Read the units parameters from the snapshot
+    units = f["InternalCodeUnits"]
+    unit_mass = units.attrs["Unit mass in cgs (U_M)"]
+    unit_length = units.attrs["Unit length in cgs (U_L)"]
+    unit_time = units.attrs["Unit time in cgs (U_t)"]
+    unit_density = unit_mass / unit_length ** 3
+
+    # Header
+    header = f["Header"]
+    BoxSize = header.attrs["BoxSize"]
+    Time = header.attrs["Time"]
+
+    # Read data
+    pos = f["/PartType0/Coordinates"][:]
+    rho = f["/PartType0/Densities"][:]
+    mass = f["/PartType0/Masses"][:]
+    ids = f["/PartType0/ParticleIDs"][:]
+
+    # Center the model and compute particle radius
+    pos = pos - BoxSize / 2
+    x = pos[:, 0]
+    y = pos[:, 1]
+    z = pos[:, 2]
+    r = np.sqrt(x * x + y * y + z * z)
+    n = len(r)
+
+    # sort particles according to their distance
+    idx = np.argsort(r)
+    r = r[idx]
+    mass = mass[idx]
+    ids = ids[idx]
+
+    nparts_per_bin = 100
+
+    # loop over bins containing nparts_per_bin particles
+    idx = np.arange(n)
+
+    # bins radius and mass
+    rs_beg = np.array([])
+    rs_end = np.array([])
+    ms = np.array([])
+
+    i = 0
+
+    while i + nparts_per_bin < n:
+
+        rs_beg = np.concatenate((rs_beg, [r[i]]))
+        rs_end = np.concatenate((rs_end, [r[i + nparts_per_bin]]))
+        m_this_bin = np.sum(mass[i : i + nparts_per_bin])
+        ms = np.concatenate((ms, [np.sum(m_this_bin)]))
+
+        # shift
+        i = i + nparts_per_bin
+
+    # compute density
+    vol = 4 / 3 * np.pi * (rs_end ** 3 - rs_beg ** 3)
+    rho = ms / vol
+
+    # compute radius, we use the mean
+    rs = 0.5 * (rs_beg + rs_end)
+
+    # convert rho to acc
+    rho = rho * unit_density / gcm3InAcc
+
+    # convert time to Myr
+    Time = Time * unit_time / MyrInSec
+
+    return rs, rho, Time
+
+
+# Do the plot
+
+plt.figure()
+
+rs, rho, Time = ComputeDensity("snapshot_0000.hdf5")
+plt.plot(rs, rho, c="b", label=r"$t=%5.1f\,\rm{[Myr]}$" % Time, lw=1)
+
+rs, rho, Time = ComputeDensity("snapshot_0050.hdf5")
+plt.plot(rs, rho, c="r", label=r"$t=%5.1f\,\rm{[Myr]}$" % Time, lw=1)
+
+plt.loglog()
+
+plt.legend()
+plt.xlabel("${\\rm{Radius~[kpc]}}$", labelpad=0)
+plt.ylabel("${\\rm{Density~[atom/cm^3]}}$", labelpad=0)
+
+
+plt.savefig("GasDensityProfile.png", dpi=200)

examples/HydroTests/NFW_Hydrostatic/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e nfw.hdf5 ]
+then
+    echo "Fetching initial conditions file for the example..."
+    ./getICs.sh
+fi
+
+
+# Run SWIFT
+../../../swift --hydro --external-gravity --self-gravity --threads=14  NFW_Hydrostatic.yml
+
+# Compute gas density profile
+python3  plotGasDensityProfile.py

examples/HydroTests/Noh_1D/makeIC.py

 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- # 
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- # GNU General Public License for more details.
- # 
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program.  If not, see <http://www.gnu.org/licenses/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
 
 import h5py
 from numpy import *
@@ -24,18 +24,18 @@ from numpy import *
 
 # Parameters
 numPart = 1001
-gamma = 5./3.      # Gas adiabatic index
-rho0 = 1.          # Background density
-P0 = 1.e-6         # Background pressure
-fileName = "noh.hdf5" 
+gamma = 5.0 / 3.0  # Gas adiabatic index
+rho0 = 1.0  # Background density
+P0 = 1.0e-6  # Background pressure
+fileName = "noh.hdf5"
 
-#---------------------------------------------------
+# ---------------------------------------------------
 coords = zeros((numPart, 3))
 h = zeros(numPart)
-vol = 2.
+vol = 2.0
 
 for i in range(numPart):
-    coords[i,0] = i * vol/numPart + vol/(2.*numPart)
+    coords[i, 0] = i * vol / numPart + vol / (2.0 * numPart)
     h[i] = 1.2348 * vol / numPart
 
 # Generate extra arrays
@@ -44,20 +44,20 @@ ids = linspace(1, numPart, numPart)
 m = zeros(numPart)
 u = zeros(numPart)
 
-m[:] = rho0 * vol / numPart    
+m[:] = rho0 * vol / numPart
 u[:] = P0 / (rho0 * (gamma - 1))
-v[coords[:,0]<vol/2. ,0] = 1
-v[coords[:,0]>vol/2. ,0] = -1
+v[coords[:, 0] < vol / 2.0, 0] = 1
+v[coords[:, 0] > vol / 2.0, 0] = -1
 
-#--------------------------------------------------
+# --------------------------------------------------
 
-#File
-file = h5py.File(fileName, 'w')
+# File
+file = h5py.File(fileName, "w")
 
 # Header
 grp = file.create_group("/Header")
 grp.attrs["BoxSize"] = [vol, vol, vol]
-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
@@ -66,21 +66,21 @@ grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 grp.attrs["Flag_Entropy_ICs"] = 0
 grp.attrs["Dimension"] = 1
 
-#Units
+# Units
 grp = file.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 1.
-grp.attrs["Unit mass in cgs (U_M)"] = 1.
-grp.attrs["Unit time in cgs (U_t)"] = 1.
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+grp.attrs["Unit length in cgs (U_L)"] = 1.0
+grp.attrs["Unit mass in cgs (U_M)"] = 1.0
+grp.attrs["Unit time in cgs (U_t)"] = 1.0
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
 
-#Particle group
+# Particle group
 grp = file.create_group("/PartType0")
-grp.create_dataset('Coordinates', data=coords, dtype='d')
-grp.create_dataset('Velocities', data=v, dtype='f')
-grp.create_dataset('Masses', data=m, dtype='f')
-grp.create_dataset('SmoothingLength', data=h, dtype='f')
-grp.create_dataset('InternalEnergy', data=u, dtype='f')
-grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+grp.create_dataset("Coordinates", data=coords, dtype="d")
+grp.create_dataset("Velocities", data=v, dtype="f")
+grp.create_dataset("Masses", data=m, dtype="f")
+grp.create_dataset("SmoothingLength", data=h, dtype="f")
+grp.create_dataset("InternalEnergy", data=u, dtype="f")
+grp.create_dataset("ParticleIDs", data=ids, dtype="L")
 
 file.close()

examples/HydroTests/Noh_1D/plotSolution.py

 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- # 
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- # GNU General Public License for more details.
- # 
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program.  If not, see <http://www.gnu.org/licenses/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
 
 # Computes the analytical solution of the Noh problem and plots the SPH answer
- 
+
 
 # Parameters
-gas_gamma = 5./3.      # Polytropic index
-rho0 = 1.              # Background density
-P0 = 1.e-6             # Background pressure
+gas_gamma = 5.0 / 3.0  # Polytropic index
+rho0 = 1.0  # Background density
+P0 = 1.0e-6  # Background pressure
 v0 = 1
 
 import matplotlib
+
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
+import numpy as np
 import h5py
+import sys
 
-# Plot parameters
-params = {'axes.labelsize': 10,
-'axes.titlesize': 10,
-'font.size': 12,
-'legend.fontsize': 12,
-'xtick.labelsize': 10,
-'ytick.labelsize': 10,
-'text.usetex': True,
- 'figure.figsize' : (9.90,6.45),
-'figure.subplot.left'    : 0.045,
-'figure.subplot.right'   : 0.99,
-'figure.subplot.bottom'  : 0.05,
-'figure.subplot.top'     : 0.99,
-'figure.subplot.wspace'  : 0.15,
-'figure.subplot.hspace'  : 0.12,
-'lines.markersize' : 6,
-'lines.linewidth' : 3.,
-'text.latex.unicode': True
-}
-rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
-
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
-
 # Read the simulation data
-sim = h5py.File("noh_%04d.hdf5"%snap, "r")
+sim = h5py.File("noh_%04d.hdf5" % snap, "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
 time = sim["/Header"].attrs["Time"][0]
 scheme = sim["/HydroScheme"].attrs["Scheme"]
@@ -67,8 +48,8 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
 eta = sim["/HydroScheme"].attrs["Kernel eta"]
 git = sim["Code"].attrs["Git Revision"]
 
-x = sim["/PartType0/Coordinates"][:,0]
-v = sim["/PartType0/Velocities"][:,0]
+x = sim["/PartType0/Coordinates"][:, 0]
+v = sim["/PartType0/Velocities"][:, 0]
 u = sim["/PartType0/InternalEnergies"][:]
 S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
@@ -83,20 +64,20 @@ x += x_min
 # ---------------------------------------------------------------
 # Don't touch anything after this.
 # ---------------------------------------------------------------
-x_s = np.arange(0, 2., 2./N) - 1.
+x_s = np.arange(0, 2.0, 2.0 / N) - 1.0
 rho_s = np.ones(N) * rho0
 P_s = np.ones(N) * rho0
 v_s = np.ones(N) * v0
 
 # Shock position
-u0 = rho0 * P0 * (gas_gamma-1)
+u0 = rho0 * P0 * (gas_gamma - 1)
 us = 0.5 * (gas_gamma - 1) * v0
 rs = us * time
 
 # Post-shock values
 rho_s[np.abs(x_s) < rs] = rho0 * (gas_gamma + 1) / (gas_gamma - 1)
-P_s[np.abs(x_s) < rs] = 0.5 * rho0 * v0**2 * (gas_gamma + 1) 
-v_s[np.abs(x_s) < rs] = 0.
+P_s[np.abs(x_s) < rs] = 0.5 * rho0 * v0 ** 2 * (gas_gamma + 1)
+v_s[np.abs(x_s) < rs] = 0.0
 
 # Pre-shock values
 rho_s[np.abs(x_s) >= rs] = rho0
@@ -105,71 +86,104 @@ v_s[x_s >= rs] = -v0
 v_s[x_s <= -rs] = v0
 
 # Additional arrays
-u_s = P_s / (rho_s * (gas_gamma - 1.))  #internal energy
-s_s = P_s / rho_s**gas_gamma # entropic function
+u_s = P_s / (rho_s * (gas_gamma - 1.0))  # internal energy
+s_s = P_s / rho_s ** gas_gamma  # entropic function
 
 # Plot the interesting quantities
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=4,
+    markeredgecolor="none",
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
+
+errorbar_props = dict(color=binned_color, ms=binned_marker_size, fmt=".", lw=1.2)
 
 # Velocity profile --------------------------------
-subplot(231)
-plot(x, v, '.', color='r', ms=4.0)
-plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_x$", labelpad=-4)
-xlim(-0.5, 0.5)
-ylim(-1.2, 1.2)
+plt.subplot(231)
+plt.plot(x, v, **scatter_props)
+plt.plot(x_s, v_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_x$", labelpad=-4)
+plt.xlim(-0.5, 0.5)
+plt.ylim(-1.2, 1.2)
 
 # Density profile --------------------------------
-subplot(232)
-plot(x, rho, '.', color='r', ms=4.0)
-plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
-xlim(-0.5, 0.5)
-ylim(0.95, 4.4)
+plt.subplot(232)
+plt.plot(x, rho, **scatter_props)
+plt.plot(x_s, rho_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0.95, 4.4)
 
 # Pressure profile --------------------------------
-subplot(233)
-plot(x, P, '.', color='r', ms=4.0)
-plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Pressure}}~P$", labelpad=0)
-xlim(-0.5, 0.5)
-ylim(-0.05, 1.8)
+plt.subplot(233)
+plt.plot(x, P, **scatter_props)
+plt.plot(x_s, P_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.xlim(-0.5, 0.5)
+plt.ylim(-0.05, 1.8)
 
 # Internal energy profile -------------------------
-subplot(234)
-plot(x, u, '.', color='r', ms=4.0)
-plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
-xlim(-0.5, 0.5)
-ylim(-0.05, 0.8)
+plt.subplot(234)
+plt.plot(x, u, **scatter_props)
+plt.plot(x_s, u_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.xlim(-0.5, 0.5)
+plt.ylim(-0.05, 0.8)
 
 # Entropy profile ---------------------------------
-subplot(235)
-plot(x, S, '.', color='r', ms=4.0)
-plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
-xlim(-0.5, 0.5)
-ylim(-0.05, 0.2)
+plt.subplot(235)
+plt.plot(x, S, **scatter_props)
+plt.plot(x_s, s_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.xlabel("${\\rm{Position}}~x$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
+plt.xlim(-0.5, 0.5)
+plt.ylim(-0.05, 0.2)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
-
-text(-0.49, 0.9, "Noh problem with  $\\gamma=%.3f$ in 1D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
-text(-0.49, 0.8, "ICs:~~ $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$"%(1e-6, 1., -1.), fontsize=10)
-plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-
-savefig("Noh.png", dpi=200)
+plt.subplot(236, frameon=False)
+
+text_fontsize = 5
+
+plt.text(
+    -0.45,
+    0.9,
+    "Noh problem with  $\\gamma=%.3f$ in 1D at $t=%.2f$" % (gas_gamma, time),
+    fontsize=text_fontsize,
+)
+plt.text(
+    -0.45,
+    0.8,
+    "ICs: $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$" % (1e-6, 1.0, -1.0),
+    fontsize=text_fontsize,
+)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
+)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("Noh.png", dpi=200)

examples/HydroTests/Noh_1D/run.sh

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

examples/HydroTests/Noh_2D/getGlass.sh

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

examples/HydroTests/Noh_2D/makeIC.py

-
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- # 
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- # GNU General Public License for more details.
- # 
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program.  If not, see <http://www.gnu.org/licenses/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
 
 import h5py
 from numpy import *
@@ -24,18 +23,18 @@ from numpy import *
 # Generates a swift IC file for the 2D Noh problem in a periodic box
 
 # Parameters
-gamma = 5./3.          # Gas adiabatic index
-gamma = 5./3.      # Gas adiabatic index
-rho0 = 1.          # Background density
-P0 = 1.e-6         # Background pressure
-fileName = "noh.hdf5" 
+gamma = 5.0 / 3.0  # Gas adiabatic index
+gamma = 5.0 / 3.0  # Gas adiabatic index
+rho0 = 1.0  # Background density
+P0 = 1.0e-6  # Background pressure
+fileName = "noh.hdf5"
 
-#---------------------------------------------------
+# ---------------------------------------------------
 glass = h5py.File("glassPlane_128.hdf5", "r")
 
-vol = 4.
+vol = 4.0
 
-pos = glass["/PartType0/Coordinates"][:,:] * sqrt(vol)
+pos = glass["/PartType0/Coordinates"][:, :] * sqrt(vol)
 h = glass["/PartType0/SmoothingLength"][:] * sqrt(vol)
 numPart = size(h)
 
@@ -45,26 +44,26 @@ ids = linspace(1, numPart, numPart)
 
 m = zeros(numPart)
 u = zeros(numPart)
-m[:] = rho0 * vol / numPart    
+m[:] = rho0 * vol / numPart
 u[:] = P0 / (rho0 * (gamma - 1))
 
 # Make radial velocities
-#r = sqrt((pos[:,0]-1.)**2 + (pos[:,1]-1.)**2)
-#theta = arctan2((pos[:,1]-1.), (pos[:,0]-1.))
-v[:,0] = -(pos[:,0] - 1.)
-v[:,1] = -(pos[:,1] - 1.)
+# r = sqrt((pos[:,0]-1.)**2 + (pos[:,1]-1.)**2)
+# theta = arctan2((pos[:,1]-1.), (pos[:,0]-1.))
+v[:, 0] = -(pos[:, 0] - 1.0)
+v[:, 1] = -(pos[:, 1] - 1.0)
 
-norm_v = sqrt(v[:,0]**2 + v[:,1]**2)
-v[:,0] /= norm_v
-v[:,1] /= norm_v
+norm_v = sqrt(v[:, 0] ** 2 + v[:, 1] ** 2)
+v[:, 0] /= norm_v
+v[:, 1] /= norm_v
 
-#File
-file = h5py.File(fileName, 'w')
+# File
+file = h5py.File(fileName, "w")
 
 # Header
 grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = [sqrt(vol), sqrt(vol), 1.]
-grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["BoxSize"] = [sqrt(vol), sqrt(vol), 1.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
@@ -73,22 +72,22 @@ grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 grp.attrs["Flag_Entropy_ICs"] = 0
 grp.attrs["Dimension"] = 2
 
-#Units
+# Units
 grp = file.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 1.
-grp.attrs["Unit mass in cgs (U_M)"] = 1.
-grp.attrs["Unit time in cgs (U_t)"] = 1.
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+grp.attrs["Unit length in cgs (U_L)"] = 1.0
+grp.attrs["Unit mass in cgs (U_M)"] = 1.0
+grp.attrs["Unit time in cgs (U_t)"] = 1.0
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
 
-#Particle group
+# Particle group
 grp = file.create_group("/PartType0")
-grp.create_dataset('Coordinates', data=pos, dtype='d')
-grp.create_dataset('Velocities', data=v, dtype='f')
-grp.create_dataset('Masses', data=m, dtype='f')
-grp.create_dataset('SmoothingLength', data=h, dtype='f')
-grp.create_dataset('InternalEnergy', data=u, dtype='f')
-grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+grp.create_dataset("Coordinates", data=pos, dtype="d")
+grp.create_dataset("Velocities", data=v, dtype="f")
+grp.create_dataset("Masses", data=m, dtype="f")
+grp.create_dataset("SmoothingLength", data=h, dtype="f")
+grp.create_dataset("InternalEnergy", data=u, dtype="f")
+grp.create_dataset("ParticleIDs", data=ids, dtype="L")
 
 
 file.close()

examples/HydroTests/Noh_2D/plotSolution.py

 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- # 
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- # GNU General Public License for more details.
- # 
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program.  If not, see <http://www.gnu.org/licenses/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
 
 # Computes the analytical solution of the Noh problem and plots the SPH answer
- 
+
 
 # Parameters
-gas_gamma = 5./3.      # Polytropic index
-rho0 = 1.              # Background density
-P0 = 1.e-6             # Background pressure
+gas_gamma = 5.0 / 3.0  # Polytropic index
+rho0 = 1.0  # Background density
+P0 = 1.0e-6  # Background pressure
 v0 = 1
 
 import matplotlib
+
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
 from scipy import stats
+import numpy as np
 import h5py
+import sys
 
-# Plot parameters
-params = {'axes.labelsize': 10,
-'axes.titlesize': 10,
-'font.size': 12,
-'legend.fontsize': 12,
-'xtick.labelsize': 10,
-'ytick.labelsize': 10,
-'text.usetex': True,
- 'figure.figsize' : (9.90,6.45),
-'figure.subplot.left'    : 0.045,
-'figure.subplot.right'   : 0.99,
-'figure.subplot.bottom'  : 0.05,
-'figure.subplot.top'     : 0.99,
-'figure.subplot.wspace'  : 0.15,
-'figure.subplot.hspace'  : 0.12,
-'lines.markersize' : 6,
-'lines.linewidth' : 3.,
-'text.latex.unicode': True
-}
-rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
-
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
 # Read the simulation data
-sim = h5py.File("noh_%04d.hdf5"%snap, "r")
+sim = h5py.File("noh_%04d.hdf5" % snap, "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
 time = sim["/Header"].attrs["Time"][0]
 scheme = sim["/HydroScheme"].attrs["Scheme"]
@@ -67,132 +49,166 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
 eta = sim["/HydroScheme"].attrs["Kernel eta"]
 git = sim["Code"].attrs["Git Revision"]
 
-x = sim["/PartType0/Coordinates"][:,0]
-y = sim["/PartType0/Coordinates"][:,1]
-vx = sim["/PartType0/Velocities"][:,0]
-vy = sim["/PartType0/Velocities"][:,1]
+x = sim["/PartType0/Coordinates"][:, 0]
+y = sim["/PartType0/Coordinates"][:, 1]
+vx = sim["/PartType0/Velocities"][:, 0]
+vy = sim["/PartType0/Velocities"][:, 1]
 u = sim["/PartType0/InternalEnergies"][:]
 S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
 rho = sim["/PartType0/Densities"][:]
 
-r = np.sqrt((x-1)**2 + (y-1)**2)
-v = -np.sqrt(vx**2 + vy**2)
-
-# Bin te data
-r_bin_edge = np.arange(0., 1., 0.02)
-r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
-rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
-v_bin,_,_ = stats.binned_statistic(r, v, statistic='mean', bins=r_bin_edge)
-P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
-S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
-u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
-rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
-v2_bin,_,_ = stats.binned_statistic(r, v**2, statistic='mean', bins=r_bin_edge)
-P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
-S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
-u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_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)
+r = np.sqrt((x - 1) ** 2 + (y - 1) ** 2)
+v = -np.sqrt(vx ** 2 + vy ** 2)
+
+# Bin the data
+r_bin_edge = np.arange(0.0, 1.0, 0.02)
+r_bin = 0.5 * (r_bin_edge[1:] + r_bin_edge[:-1])
+rho_bin, _, _ = stats.binned_statistic(r, rho, statistic="mean", bins=r_bin_edge)
+v_bin, _, _ = stats.binned_statistic(r, v, statistic="mean", bins=r_bin_edge)
+P_bin, _, _ = stats.binned_statistic(r, P, statistic="mean", bins=r_bin_edge)
+S_bin, _, _ = stats.binned_statistic(r, S, statistic="mean", bins=r_bin_edge)
+u_bin, _, _ = stats.binned_statistic(r, u, statistic="mean", bins=r_bin_edge)
+rho2_bin, _, _ = stats.binned_statistic(r, rho ** 2, statistic="mean", bins=r_bin_edge)
+v2_bin, _, _ = stats.binned_statistic(r, v ** 2, statistic="mean", bins=r_bin_edge)
+P2_bin, _, _ = stats.binned_statistic(r, P ** 2, statistic="mean", bins=r_bin_edge)
+S2_bin, _, _ = stats.binned_statistic(r, S ** 2, statistic="mean", bins=r_bin_edge)
+u2_bin, _, _ = stats.binned_statistic(r, u ** 2, statistic="mean", bins=r_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)
 
 # Analytic solution
 N = 1000  # Number of points
 
-x_s = np.arange(0, 2., 2./N) - 1.
+x_s = np.arange(0, 2.0, 2.0 / N) - 1.0
 rho_s = np.ones(N) * rho0
 P_s = np.ones(N) * rho0
 v_s = np.ones(N) * v0
 
 # Shock position
-u0 = rho0 * P0 * (gas_gamma-1)
+u0 = rho0 * P0 * (gas_gamma - 1)
 us = 0.5 * (gas_gamma - 1) * v0
 rs = us * time
 
 # Post-shock values
-rho_s[np.abs(x_s) < rs] = rho0 * ((gas_gamma + 1) / (gas_gamma - 1))**2
-P_s[np.abs(x_s) < rs] = 0.5 * rho0 * v0**2 * (gas_gamma + 1)**2 / (gas_gamma-1)
-v_s[np.abs(x_s) < rs] = 0.
+rho_s[np.abs(x_s) < rs] = rho0 * ((gas_gamma + 1) / (gas_gamma - 1)) ** 2
+P_s[np.abs(x_s) < rs] = 0.5 * rho0 * v0 ** 2 * (gas_gamma + 1) ** 2 / (gas_gamma - 1)
+v_s[np.abs(x_s) < rs] = 0.0
 
 # Pre-shock values
-rho_s[np.abs(x_s) >= rs] = rho0 * (1 + v0 * time/np.abs(x_s[np.abs(x_s) >=rs]))
+rho_s[np.abs(x_s) >= rs] = rho0 * (1 + v0 * time / np.abs(x_s[np.abs(x_s) >= rs]))
 P_s[np.abs(x_s) >= rs] = 0
 v_s[x_s >= rs] = -v0
 v_s[x_s <= -rs] = v0
 
 # Additional arrays
-u_s = P_s / (rho_s * (gas_gamma - 1.))  #internal energy
-s_s = P_s / rho_s**gas_gamma # entropic function
+u_s = P_s / (rho_s * (gas_gamma - 1.0))  # internal energy
+s_s = P_s / rho_s ** gas_gamma  # entropic function
 
 # Plot the interesting quantities
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=1,
+    markeredgecolor="none",
+    alpha=0.2,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
+
+errorbar_props = dict(color=binned_color, ms=binned_marker_size, fmt=".", lw=1.2)
 
 # Velocity profile --------------------------------
-subplot(231)
-plot(r, v, '.', color='r', ms=0.5, alpha=0.2)
-plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(r_bin, v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_r$", labelpad=-4)
-xlim(0, 0.5)
-ylim(-1.2, 0.4)
+plt.subplot(231)
+plt.plot(r, v, **scatter_props)
+plt.plot(x_s, v_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, v_bin, yerr=v_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_r$", labelpad=-4)
+plt.xlim(0, 0.5)
+plt.ylim(-1.2, 0.4)
 
 # Density profile --------------------------------
-subplot(232)
-plot(r, rho, '.', color='r', ms=0.5, alpha=0.2)
-plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
-xlim(0, 0.5)
-ylim(0.95, 19)
+plt.subplot(232)
+plt.plot(r, rho, **scatter_props)
+plt.plot(x_s, rho_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.xlim(0, 0.5)
+plt.ylim(0.95, 19)
 
 # Pressure profile --------------------------------
-subplot(233)
-plot(r, P, '.', color='r', ms=0.5, alpha=0.2)
-plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Pressure}}~P$", labelpad=0)
-xlim(0, 0.5)
-ylim(-0.5, 11)
+plt.subplot(233)
+plt.plot(r, P, **scatter_props)
+plt.plot(x_s, P_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, P_bin, yerr=P_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.xlim(0, 0.5)
+plt.ylim(-0.5, 11)
 
 # Internal energy profile -------------------------
-subplot(234)
-plot(r, u, '.', color='r', ms=0.5, alpha=0.2)
-plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
-xlim(0, 0.5)
-ylim(-0.05, 0.8)
+plt.subplot(234)
+plt.plot(r, u, **scatter_props)
+plt.plot(x_s, u_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, u_bin, yerr=u_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.xlim(0, 0.5)
+plt.ylim(-0.05, 0.8)
 
 # Entropy profile ---------------------------------
-subplot(235)
-plot(r, S, '.', color='r', ms=0.5, alpha=0.2) 
-plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
-xlim(0, 0.5)
-ylim(-0.05, 0.2)
+plt.subplot(235)
+plt.plot(r, S, **scatter_props)
+plt.plot(x_s, s_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, S_bin, yerr=S_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
+plt.xlim(0, 0.5)
+plt.ylim(-0.05, 0.2)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
-
-text(-0.49, 0.9, "Noh problem with  $\\gamma=%.3f$ in 2D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
-text(-0.49, 0.8, "ICs:~~ $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$"%(1e-6, 1., -1.), fontsize=10)
-plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-
-savefig("Noh.png", dpi=200)
+plt.subplot(236, frameon=False)
+
+text_fontsize = 5
+
+plt.text(
+    -0.45,
+    0.9,
+    "Noh problem with  $\\gamma=%.3f$ in 2D at $t=%.2f$" % (gas_gamma, time),
+    fontsize=text_fontsize,
+)
+plt.text(
+    -0.45,
+    0.8,
+    "ICs: $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$" % (1e-6, 1.0, -1.0),
+    fontsize=text_fontsize,
+)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
+)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("Noh.png", dpi=200)

examples/HydroTests/Noh_2D/run.sh

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

examples/HydroTests/Noh_3D/getGlass.sh

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

examples/HydroTests/Noh_3D/makeIC.py

-
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- # 
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- # GNU General Public License for more details.
- # 
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program.  If not, see <http://www.gnu.org/licenses/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
 
 import h5py
 from numpy import *
@@ -24,19 +23,19 @@ from numpy import *
 # Generates a swift IC file for the 3D Noh problem in a periodic box
 
 # Parameters
-gamma = 5./3.          # Gas adiabatic index
-gamma = 5./3.      # Gas adiabatic index
-rho0 = 1.          # Background density
-P0 = 1.e-6         # Background pressure
-fileName = "noh.hdf5" 
+gamma = 5.0 / 3.0  # Gas adiabatic index
+gamma = 5.0 / 3.0  # Gas adiabatic index
+rho0 = 1.0  # Background density
+P0 = 1.0e-6  # Background pressure
+fileName = "noh.hdf5"
 
-#---------------------------------------------------
+# ---------------------------------------------------
 glass = h5py.File("glassCube_64.hdf5", "r")
 
-vol = 8.
+vol = 8.0
 
-pos = glass["/PartType0/Coordinates"][:,:] * vol**(1./3.)
-h = glass["/PartType0/SmoothingLength"][:] * vol**(1./3.)
+pos = glass["/PartType0/Coordinates"][:, :] * vol ** (1.0 / 3.0)
+h = glass["/PartType0/SmoothingLength"][:] * vol ** (1.0 / 3.0)
 numPart = size(h)
 
 # Generate extra arrays
@@ -45,28 +44,28 @@ ids = linspace(1, numPart, numPart)
 
 m = zeros(numPart)
 u = zeros(numPart)
-m[:] = rho0 * vol / numPart    
+m[:] = rho0 * vol / numPart
 u[:] = P0 / (rho0 * (gamma - 1))
 
 # Make radial velocities
-#r = sqrt((pos[:,0]-1.)**2 + (pos[:,1]-1.)**2)
-#theta = arctan2((pos[:,1]-1.), (pos[:,0]-1.))
-v[:,0] = -(pos[:,0] - 1.)
-v[:,1] = -(pos[:,1] - 1.)
-v[:,2] = -(pos[:,2] - 1.)
+# r = sqrt((pos[:,0]-1.)**2 + (pos[:,1]-1.)**2)
+# theta = arctan2((pos[:,1]-1.), (pos[:,0]-1.))
+v[:, 0] = -(pos[:, 0] - 1.0)
+v[:, 1] = -(pos[:, 1] - 1.0)
+v[:, 2] = -(pos[:, 2] - 1.0)
 
-norm_v = sqrt(v[:,0]**2 + v[:,1]**2 + v[:,2]**2)
-v[:,0] /= norm_v
-v[:,1] /= norm_v
-v[:,2] /= norm_v
+norm_v = sqrt(v[:, 0] ** 2 + v[:, 1] ** 2 + v[:, 2] ** 2)
+v[:, 0] /= norm_v
+v[:, 1] /= norm_v
+v[:, 2] /= norm_v
 
-#File
-file = h5py.File(fileName, 'w')
+# File
+file = h5py.File(fileName, "w")
 
 # Header
 grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = [vol**(1./3.), vol**(1./3.), vol**(1./3.)]
-grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["BoxSize"] = [vol ** (1.0 / 3.0), vol ** (1.0 / 3.0), vol ** (1.0 / 3.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
@@ -75,22 +74,22 @@ grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 grp.attrs["Flag_Entropy_ICs"] = 0
 grp.attrs["Dimension"] = 3
 
-#Units
+# Units
 grp = file.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 1.
-grp.attrs["Unit mass in cgs (U_M)"] = 1.
-grp.attrs["Unit time in cgs (U_t)"] = 1.
-grp.attrs["Unit current in cgs (U_I)"] = 1.
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+grp.attrs["Unit length in cgs (U_L)"] = 1.0
+grp.attrs["Unit mass in cgs (U_M)"] = 1.0
+grp.attrs["Unit time in cgs (U_t)"] = 1.0
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
 
-#Particle group
+# Particle group
 grp = file.create_group("/PartType0")
-grp.create_dataset('Coordinates', data=pos, dtype='d')
-grp.create_dataset('Velocities', data=v, dtype='f')
-grp.create_dataset('Masses', data=m, dtype='f')
-grp.create_dataset('SmoothingLength', data=h, dtype='f')
-grp.create_dataset('InternalEnergy', data=u, dtype='f')
-grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+grp.create_dataset("Coordinates", data=pos, dtype="d")
+grp.create_dataset("Velocities", data=v, dtype="f")
+grp.create_dataset("Masses", data=m, dtype="f")
+grp.create_dataset("SmoothingLength", data=h, dtype="f")
+grp.create_dataset("InternalEnergy", data=u, dtype="f")
+grp.create_dataset("ParticleIDs", data=ids, dtype="L")
 
 
 file.close()

examples/HydroTests/Noh_3D/plotSolution.py

 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- # 
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- # GNU General Public License for more details.
- # 
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program.  If not, see <http://www.gnu.org/licenses/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
 
 # Computes the analytical solution of the Noh problem and plots the SPH answer
- 
+
 
 # Parameters
-gas_gamma = 5./3.      # Polytropic index
-rho0 = 1.              # Background density
-P0 = 1.e-6             # Background pressure
+gas_gamma = 5.0 / 3.0  # Polytropic index
+rho0 = 1.0  # Background density
+P0 = 1.0e-6  # Background pressure
 v0 = 1
 
 import matplotlib
+
 matplotlib.use("Agg")
-from pylab import *
+import matplotlib.pyplot as plt
 from scipy import stats
+import numpy as np
 import h5py
+import sys
 
-# Plot parameters
-params = {'axes.labelsize': 10,
-'axes.titlesize': 10,
-'font.size': 12,
-'legend.fontsize': 12,
-'xtick.labelsize': 10,
-'ytick.labelsize': 10,
-'text.usetex': True,
- 'figure.figsize' : (9.90,6.45),
-'figure.subplot.left'    : 0.045,
-'figure.subplot.right'   : 0.99,
-'figure.subplot.bottom'  : 0.05,
-'figure.subplot.top'     : 0.99,
-'figure.subplot.wspace'  : 0.15,
-'figure.subplot.hspace'  : 0.12,
-'lines.markersize' : 6,
-'lines.linewidth' : 3.,
-'text.latex.unicode': True
-}
-rcParams.update(params)
-rc('font',**{'family':'sans-serif','sans-serif':['Times']})
-
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
 
 snap = int(sys.argv[1])
 
-
 # Read the simulation data
-sim = h5py.File("noh_%04d.hdf5"%snap, "r")
+sim = h5py.File("noh_%04d.hdf5" % snap, "r")
 boxSize = sim["/Header"].attrs["BoxSize"][0]
 time = sim["/Header"].attrs["Time"][0]
 scheme = sim["/HydroScheme"].attrs["Scheme"]
@@ -68,135 +49,171 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
 eta = sim["/HydroScheme"].attrs["Kernel eta"]
 git = sim["Code"].attrs["Git Revision"]
 
-x = sim["/PartType0/Coordinates"][:,0]
-y = sim["/PartType0/Coordinates"][:,1]
-z = sim["/PartType0/Coordinates"][:,2]
-vx = sim["/PartType0/Velocities"][:,0]
-vy = sim["/PartType0/Velocities"][:,1]
-vz = sim["/PartType0/Velocities"][:,2]
+x = sim["/PartType0/Coordinates"][:, 0]
+y = sim["/PartType0/Coordinates"][:, 1]
+z = sim["/PartType0/Coordinates"][:, 2]
+vx = sim["/PartType0/Velocities"][:, 0]
+vy = sim["/PartType0/Velocities"][:, 1]
+vz = sim["/PartType0/Velocities"][:, 2]
 u = sim["/PartType0/InternalEnergies"][:]
 S = sim["/PartType0/Entropies"][:]
 P = sim["/PartType0/Pressures"][:]
 rho = sim["/PartType0/Densities"][:]
 
-r = np.sqrt((x-1)**2 + (y-1)**2 + (z-1)**2)
-v = -np.sqrt(vx**2 + vy**2 + vz**2)
-
-# Bin te data
-r_bin_edge = np.arange(0., 1., 0.02)
-r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
-rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
-v_bin,_,_ = stats.binned_statistic(r, v, statistic='mean', bins=r_bin_edge)
-P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
-S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
-u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
-rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
-v2_bin,_,_ = stats.binned_statistic(r, v**2, statistic='mean', bins=r_bin_edge)
-P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
-S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
-u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_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)
+r = np.sqrt((x - 1) ** 2 + (y - 1) ** 2 + (z - 1) ** 2)
+v = -np.sqrt(vx ** 2 + vy ** 2 + vz ** 2)
+
+# Bin the data
+r_bin_edge = np.arange(0.0, 1.0, 0.02)
+r_bin = 0.5 * (r_bin_edge[1:] + r_bin_edge[:-1])
+rho_bin, _, _ = stats.binned_statistic(r, rho, statistic="mean", bins=r_bin_edge)
+v_bin, _, _ = stats.binned_statistic(r, v, statistic="mean", bins=r_bin_edge)
+P_bin, _, _ = stats.binned_statistic(r, P, statistic="mean", bins=r_bin_edge)
+S_bin, _, _ = stats.binned_statistic(r, S, statistic="mean", bins=r_bin_edge)
+u_bin, _, _ = stats.binned_statistic(r, u, statistic="mean", bins=r_bin_edge)
+rho2_bin, _, _ = stats.binned_statistic(r, rho ** 2, statistic="mean", bins=r_bin_edge)
+v2_bin, _, _ = stats.binned_statistic(r, v ** 2, statistic="mean", bins=r_bin_edge)
+P2_bin, _, _ = stats.binned_statistic(r, P ** 2, statistic="mean", bins=r_bin_edge)
+S2_bin, _, _ = stats.binned_statistic(r, S ** 2, statistic="mean", bins=r_bin_edge)
+u2_bin, _, _ = stats.binned_statistic(r, u ** 2, statistic="mean", bins=r_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)
 
 
 # Analytic solution
 N = 1000  # Number of points
 
-x_s = np.arange(0, 2., 2./N) - 1.
+x_s = np.arange(0, 2.0, 2.0 / N) - 1.0
 rho_s = np.ones(N) * rho0
 P_s = np.ones(N) * rho0
 v_s = np.ones(N) * v0
 
 # Shock position
-u0 = rho0 * P0 * (gas_gamma-1)
+u0 = rho0 * P0 * (gas_gamma - 1)
 us = 0.5 * (gas_gamma - 1) * v0
 rs = us * time
 
 # Post-shock values
-rho_s[np.abs(x_s) < rs] = rho0 * ((gas_gamma + 1) / (gas_gamma - 1))**3
-P_s[np.abs(x_s) < rs] = 0.5 * rho0 * v0**2 * (gas_gamma + 1)**3 / (gas_gamma-1)**2
-v_s[np.abs(x_s) < rs] = 0.
+rho_s[np.abs(x_s) < rs] = rho0 * ((gas_gamma + 1) / (gas_gamma - 1)) ** 3
+P_s[np.abs(x_s) < rs] = (
+    0.5 * rho0 * v0 ** 2 * (gas_gamma + 1) ** 3 / (gas_gamma - 1) ** 2
+)
+v_s[np.abs(x_s) < rs] = 0.0
 
 # Pre-shock values
-rho_s[np.abs(x_s) >= rs] = rho0 * (1 + v0 * time/np.abs(x_s[np.abs(x_s) >=rs]))**2
+rho_s[np.abs(x_s) >= rs] = rho0 * (1 + v0 * time / np.abs(x_s[np.abs(x_s) >= rs])) ** 2
 P_s[np.abs(x_s) >= rs] = 0
 v_s[x_s >= rs] = -v0
 v_s[x_s <= -rs] = v0
 
 # Additional arrays
-u_s = P_s / (rho_s * (gas_gamma - 1.))  #internal energy
-s_s = P_s / rho_s**gas_gamma # entropic function
+u_s = P_s / (rho_s * (gas_gamma - 1.0))  # internal energy
+s_s = P_s / rho_s ** gas_gamma  # entropic function
 
 # Plot the interesting quantities
-figure()
+plt.figure(figsize=(7, 7 / 1.6))
+
+line_color = "C4"
+binned_color = "C2"
+binned_marker_size = 4
+
+scatter_props = dict(
+    marker=".",
+    ms=1,
+    markeredgecolor="none",
+    alpha=0.2,
+    zorder=-1,
+    rasterized=True,
+    linestyle="none",
+)
+
+errorbar_props = dict(color=binned_color, ms=binned_marker_size, fmt=".", lw=1.2)
 
 # Velocity profile --------------------------------
-subplot(231)
-plot(r, v, '.', color='r', ms=0.5, alpha=0.2)
-plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(r_bin, v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_r$", labelpad=-4)
-xlim(0, 0.5)
-ylim(-1.2, 0.4)
+plt.subplot(231)
+plt.plot(r, v, **scatter_props)
+plt.plot(x_s, v_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, v_bin, yerr=v_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Velocity}}~v_r$", labelpad=-4)
+plt.xlim(0, 0.5)
+plt.ylim(-1.2, 0.4)
 
 # Density profile --------------------------------
-subplot(232)
-plot(r, rho, '.', color='r', ms=0.5, alpha=0.2)
-plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
-xlim(0, 0.5)
-ylim(0.95, 71)
+plt.subplot(232)
+plt.plot(r, rho, **scatter_props)
+plt.plot(x_s, rho_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+plt.xlim(0, 0.5)
+plt.ylim(0.95, 71)
 
 # Pressure profile --------------------------------
-subplot(233)
-plot(r, P, '.', color='r', ms=0.5, alpha=0.2)
-plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Pressure}}~P$", labelpad=0)
-xlim(0, 0.5)
-ylim(-0.5, 25)
+plt.subplot(233)
+plt.plot(r, P, **scatter_props)
+plt.plot(x_s, P_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, P_bin, yerr=P_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+plt.xlim(0, 0.5)
+plt.ylim(-0.5, 25)
 
 # Internal energy profile -------------------------
-subplot(234)
-plot(r, u, '.', color='r', ms=0.5, alpha=0.2)
-plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
-xlim(0, 0.5)
-ylim(-0.05, 0.8)
+plt.subplot(234)
+plt.plot(r, u, **scatter_props)
+plt.plot(x_s, u_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, u_bin, yerr=u_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+plt.xlim(0, 0.5)
+plt.ylim(-0.05, 0.8)
 
 # Entropy profile ---------------------------------
-subplot(235)
-plot(r, S, '.', color='r', ms=0.5, alpha=0.2) 
-plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
-errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
-xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
-xlim(0, 0.5)
-ylim(-0.05, 0.2)
+plt.subplot(235)
+plt.plot(r, S, **scatter_props)
+plt.plot(x_s, s_s, "--", color=line_color, alpha=0.8, lw=1.2)
+plt.errorbar(r_bin, S_bin, yerr=S_sigma_bin, **errorbar_props)
+plt.xlabel("${\\rm{Radius}}~r$", labelpad=0)
+plt.ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
+plt.xlim(0, 0.5)
+plt.ylim(-0.05, 0.2)
 
 # Information -------------------------------------
-subplot(236, frameon=False)
-
-text(-0.49, 0.9, "Noh problem with  $\\gamma=%.3f$ in 3D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
-text(-0.49, 0.8, "ICs:~~ $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$"%(1e-6, 1., -1.), fontsize=10)
-plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
-text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
-text(-0.49, 0.4, scheme, fontsize=10)
-text(-0.49, 0.3, kernel, fontsize=10)
-text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
-xlim(-0.5, 0.5)
-ylim(0, 1)
-xticks([])
-yticks([])
-
-
-savefig("Noh.png", dpi=200)
+plt.subplot(236, frameon=False)
+
+text_fontsize = 5
+
+plt.text(
+    -0.45,
+    0.9,
+    "Noh problem with  $\\gamma=%.3f$ in 3D at $t=%.2f$" % (gas_gamma, time),
+    fontsize=text_fontsize,
+)
+plt.text(
+    -0.45,
+    0.8,
+    "ICs: $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$" % (1e-6, 1.0, -1.0),
+    fontsize=text_fontsize,
+)
+plt.plot([-0.45, 0.1], [0.62, 0.62], "k-", lw=1)
+plt.text(-0.45, 0.5, "$SWIFT$ %s" % git.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.4, scheme.decode("utf-8"), fontsize=text_fontsize)
+plt.text(-0.45, 0.3, kernel.decode("utf-8"), fontsize=text_fontsize)
+plt.text(
+    -0.45,
+    0.2,
+    "$%.2f$ neighbours ($\\eta=%.3f$)" % (neighbours, eta),
+    fontsize=text_fontsize,
+)
+plt.xlim(-0.5, 0.5)
+plt.ylim(0, 1)
+plt.xticks([])
+plt.yticks([])
+
+plt.tight_layout()
+
+plt.savefig("Noh.png", dpi=200)

examples/HydroTests/Noh_3D/run.sh

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