diff --git a/examples/GravityTests/HydrostaticHalo/makeIC.py b/examples/GravityTests/HydrostaticHalo/makeIC.py
index b8a4036b77c430866f700047fd06bf2c8de490e7..4b9b1488765d4c4f9db49048954e541fbb49298a 100644
--- a/examples/GravityTests/HydrostaticHalo/makeIC.py
+++ b/examples/GravityTests/HydrostaticHalo/makeIC.py
@@ -62,15 +62,15 @@ M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
const_unit_mass_in_cgs = M_vir_cgs
const_unit_length_in_cgs = r_vir_cgs
-print "UnitMass_in_cgs: ", const_unit_mass_in_cgs
-print "UnitLength_in_cgs: ", const_unit_length_in_cgs
-print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
+print("UnitMass_in_cgs: ", const_unit_mass_in_cgs)
+print("UnitLength_in_cgs: ", const_unit_length_in_cgs)
+print("UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs)
#derived quantities
const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
-print "UnitTime_in_cgs: ", const_unit_time_in_cgs
+print("UnitTime_in_cgs: ", const_unit_time_in_cgs)
const_G = ((CONST_G_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
-print 'G=', const_G
+print('G=', const_G)
# Parameters
periodic= 1 # 1 For periodic box
@@ -108,9 +108,9 @@ coords[:,2] = radius * ctheta
#shift to centre of box
coords += np.full((N,3),boxSize/2.)
-print "x range = (%f,%f)" %(np.min(coords[:,0]),np.max(coords[:,0]))
-print "y range = (%f,%f)" %(np.min(coords[:,1]),np.max(coords[:,1]))
-print "z range = (%f,%f)" %(np.min(coords[:,2]),np.max(coords[:,2]))
+print("x range = (%f,%f)" %(np.min(coords[:,0]),np.max(coords[:,0])))
+print("y range = (%f,%f)" %(np.min(coords[:,1]),np.max(coords[:,1])))
+print("z range = (%f,%f)" %(np.min(coords[:,2]),np.max(coords[:,2])))
#print np.mean(coords[:,0])
#print np.mean(coords[:,1])
@@ -156,7 +156,7 @@ z_coords = z_coords[ind]
N = x_coords.size
-print "Number of particles in the box = " , N
+print("Number of particles in the box = " , N)
#make the coords and radius arrays again
coords_2 = np.zeros((N,3))
@@ -168,13 +168,13 @@ radius = np.sqrt(coords_2[:,0]**2 + coords_2[:,1]**2 + coords_2[:,2]**2)
#test we've done it right
-print "x range = (%f,%f)" %(np.min(coords_2[:,0]),np.max(coords_2[:,0]))
-print "y range = (%f,%f)" %(np.min(coords_2[:,1]),np.max(coords_2[:,1]))
-print "z range = (%f,%f)" %(np.min(coords_2[:,2]),np.max(coords_2[:,2]))
+print("x range = (%f,%f)" %(np.min(coords_2[:,0]),np.max(coords_2[:,0])))
+print("y range = (%f,%f)" %(np.min(coords_2[:,1]),np.max(coords_2[:,1])))
+print("z range = (%f,%f)" %(np.min(coords_2[:,2]),np.max(coords_2[:,2])))
-print np.mean(coords_2[:,0])
-print np.mean(coords_2[:,1])
-print np.mean(coords_2[:,2])
+print(np.mean(coords_2[:,0]))
+print(np.mean(coords_2[:,1]))
+print(np.mean(coords_2[:,2]))
# Header
grp = file.create_group("/Header")
diff --git a/examples/HydroTests/SodShock_BCC_3D/README b/examples/HydroTests/SodShock_BCC_3D/README
new file mode 100644
index 0000000000000000000000000000000000000000..b836f33c6024a0de626d24b71309da2e5cccdbcd
--- /dev/null
+++ b/examples/HydroTests/SodShock_BCC_3D/README
@@ -0,0 +1,7 @@
+SodShock BCC 3D
+---------------
+
+This example is the exact same as the SodShock_3D, but instead
+of using a glass file, it uses a 'fake glass' that is generated
+as an exact body centered cubic lattice with the makeGlass.py
+script.
diff --git a/examples/HydroTests/SodShock_BCC_3D/analyticSolution.py b/examples/HydroTests/SodShock_BCC_3D/analyticSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..2cac6b0fabde207c325a1969434b6de9e90df9ca
--- /dev/null
+++ b/examples/HydroTests/SodShock_BCC_3D/analyticSolution.py
@@ -0,0 +1,179 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+# 2019 Josh Borrow (joshua.boorrow@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/>.
+#
+##############################################################################
+
+from numpy import *
+
+
+def analytic(
+ time, # Simulation t
+ gas_gamma=5.0 / 3.0, # Polytropic index
+ rho_L=1.0, # Density left state
+ rho_R=0.125, # Density right state
+ v_L=0.0, # Velocity left state
+ v_R=0.0, # Velocity right state
+ P_L=1.0, # Pressure left state
+ P_R=0.1, # Pressure right state
+ N=1000,
+ x_min=-1.0,
+ x_max=1.0,
+):
+ # Analytic solution
+ c_L = sqrt(gas_gamma * P_L / rho_L) # Speed of the rarefaction wave
+ c_R = sqrt(gas_gamma * P_R / rho_R) # Speed of the shock front
+
+ # Helpful variable
+ Gama = (gas_gamma - 1.0) / (gas_gamma + 1.0)
+ beta = (gas_gamma - 1.0) / (2.0 * gas_gamma)
+
+ # Characteristic function and its derivative, following Toro (2009)
+ def compute_f(P_3, P, c):
+ u = P_3 / P
+ if u > 1:
+ term1 = gas_gamma * ((gas_gamma + 1.0) * u + gas_gamma - 1.0)
+ term2 = sqrt(2.0 / term1)
+ fp = (u - 1.0) * c * term2
+ dfdp = (
+ c * term2 / P
+ + (u - 1.0)
+ * c
+ / term2
+ * (-1.0 / term1 ** 2)
+ * gas_gamma
+ * (gas_gamma + 1.0)
+ / P
+ )
+ else:
+ fp = (u ** beta - 1.0) * (2.0 * c / (gas_gamma - 1.0))
+ dfdp = 2.0 * c / (gas_gamma - 1.0) * beta * u ** (beta - 1.0) / P
+ return (fp, dfdp)
+
+ # Solution of the Riemann problem following Toro (2009)
+ def RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R):
+ P_new = (
+ (c_L + c_R + (v_L - v_R) * 0.5 * (gas_gamma - 1.0))
+ / (c_L / P_L ** beta + c_R / P_R ** beta)
+ ) ** (1.0 / beta)
+ P_3 = 0.5 * (P_R + P_L)
+ f_L = 1.0
+ while fabs(P_3 - P_new) > 1e-6:
+ P_3 = P_new
+ (f_L, dfdp_L) = compute_f(P_3, P_L, c_L)
+ (f_R, dfdp_R) = compute_f(P_3, P_R, c_R)
+ f = f_L + f_R + (v_R - v_L)
+ df = dfdp_L + dfdp_R
+ dp = -f / df
+ prnew = P_3 + dp
+ v_3 = v_L - f_L
+ return (P_new, v_3)
+
+ # Solve Riemann problem for post-shock region
+ (P_3, v_3) = RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R)
+
+ # Check direction of shocks and wave
+ shock_R = P_3 > P_R
+ shock_L = P_3 > P_L
+
+ # Velocity of shock front and and rarefaction wave
+ if shock_R:
+ v_right = v_R + c_R ** 2 * (P_3 / P_R - 1.0) / (gas_gamma * (v_3 - v_R))
+ else:
+ v_right = c_R + 0.5 * (gas_gamma + 1.0) * v_3 - 0.5 * (gas_gamma - 1.0) * v_R
+
+ if shock_L:
+ v_left = v_L + c_L ** 2 * (P_3 / p_L - 1.0) / (gas_gamma * (v_3 - v_L))
+ else:
+ v_left = c_L - 0.5 * (gas_gamma + 1.0) * v_3 + 0.5 * (gas_gamma - 1.0) * v_L
+
+ # Compute position of the transitions
+ x_23 = -fabs(v_left) * time
+ if shock_L:
+ x_12 = -fabs(v_left) * time
+ else:
+ x_12 = -(c_L - v_L) * time
+
+ x_34 = v_3 * time
+
+ x_45 = fabs(v_right) * time
+ if shock_R:
+ x_56 = fabs(v_right) * time
+ else:
+ x_56 = (c_R + v_R) * time
+
+ # Prepare arrays
+ delta_x = (x_max - x_min) / N
+ x_s = arange(x_min, x_max, delta_x)
+ rho_s = zeros(N)
+ P_s = zeros(N)
+ v_s = zeros(N)
+
+ # Compute solution in the different regions
+ for i in range(N):
+ if x_s[i] <= x_12:
+ rho_s[i] = rho_L
+ P_s[i] = P_L
+ v_s[i] = v_L
+ if x_s[i] >= x_12 and x_s[i] < x_23:
+ if shock_L:
+ rho_s[i] = rho_L * (Gama + P_3 / P_L) / (1.0 + Gama * P_3 / P_L)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ else:
+ rho_s[i] = rho_L * (
+ Gama * (0.0 - x_s[i]) / (c_L * time)
+ + Gama * v_L / c_L
+ + (1.0 - Gama)
+ ) ** (2.0 / (gas_gamma - 1.0))
+ P_s[i] = P_L * (rho_s[i] / rho_L) ** gas_gamma
+ v_s[i] = (1.0 - Gama) * (c_L - (0.0 - x_s[i]) / time) + Gama * v_L
+ if x_s[i] >= x_23 and x_s[i] < x_34:
+ if shock_L:
+ rho_s[i] = rho_L * (Gama + P_3 / P_L) / (1 + Gama * P_3 / p_L)
+ else:
+ rho_s[i] = rho_L * (P_3 / P_L) ** (1.0 / gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_34 and x_s[i] < x_45:
+ if shock_R:
+ rho_s[i] = rho_R * (Gama + P_3 / P_R) / (1.0 + Gama * P_3 / P_R)
+ else:
+ rho_s[i] = rho_R * (P_3 / P_R) ** (1.0 / gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_45 and x_s[i] < x_56:
+ if shock_R:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+ else:
+ rho_s[i] = rho_R * (
+ Gama * (x_s[i]) / (c_R * time) - Gama * v_R / c_R + (1.0 - Gama)
+ ) ** (2.0 / (gas_gamma - 1.0))
+ P_s[i] = p_R * (rho_s[i] / rho_R) ** gas_gamma
+ v_s[i] = (1.0 - Gama) * (-c_R - (-x_s[i]) / time) + Gama * v_R
+ if x_s[i] >= x_56:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+
+ # Additional arrays
+ u_s = P_s / (rho_s * (gas_gamma - 1.0)) # internal energy
+ s_s = P_s / rho_s ** gas_gamma # entropic function
+
+ return dict(x=x_s + 1, v=v_s, rho=rho_s, P=P_s, u=u_s, S=s_s)
diff --git a/examples/HydroTests/SodShock_BCC_3D/makeGlass.py b/examples/HydroTests/SodShock_BCC_3D/makeGlass.py
new file mode 100644
index 0000000000000000000000000000000000000000..974c445668e4f6a93bd6328e2f5f2c953e8ef1a6
--- /dev/null
+++ b/examples/HydroTests/SodShock_BCC_3D/makeGlass.py
@@ -0,0 +1,90 @@
+"""
+Create a fcc glass.
+"""
+
+import numpy as np
+import h5py
+
+from unyt import cm, g, s, erg
+from unyt.unit_systems import cgs_unit_system
+
+from swiftsimio import Writer
+
+
+def generate_cube(num_on_side, side_length=1.0):
+ """
+ Generates a cube
+ """
+
+ values = np.linspace(0.0, side_length, num_on_side + 1)[:-1]
+
+ positions = np.empty((num_on_side**3, 3), dtype=float)
+
+ for x in range(num_on_side):
+ for y in range(num_on_side):
+ for z in range(num_on_side):
+ index = x * num_on_side + y * num_on_side**2 + z
+
+ positions[index, 0] = values[x]
+ positions[index, 1] = values[y]
+ positions[index, 2] = values[z]
+
+ return positions
+
+
+def generate_two_cube(num_on_side, side_length=1.0):
+ cube = generate_cube(num_on_side // 2, side_length)
+
+ mips = side_length / num_on_side
+
+ positions = np.concatenate([cube, cube + mips * 0.5])
+
+ return positions
+
+
+def write_out_glass(filename, cube, side_length=1.0):
+ x = Writer(cgs_unit_system, side_length * cm)
+
+ x.gas.coordinates = cube * cm
+
+ x.gas.velocities = np.zeros_like(cube) * cm / s
+
+ x.gas.masses = np.ones(cube.shape[0], dtype=float) * g
+
+ x.gas.internal_energy = np.ones(cube.shape[0], dtype=float) * erg / g
+
+ x.gas.generate_smoothing_lengths(boxsize=side_length * cm, dimension=3)
+
+ x.write(filename)
+
+ return
+
+
+if __name__ == "__main__":
+ import argparse as ap
+
+ parser = ap.ArgumentParser(
+ description="Generate a glass file with a FCC lattice"
+ )
+
+ parser.add_argument(
+ "-n",
+ "--numparts",
+ help="Number of particles on a side. Default: 64",
+ default=1,
+ type=int
+ )
+
+ parser.add_argument(
+ "-o",
+ "--output",
+ help="Output filename.",
+ type=str,
+ required=False
+ )
+
+ args = parser.parse_args()
+
+ glass_cube = generate_two_cube(args.numparts)
+ write_out_glass(args.output, glass_cube)
+
diff --git a/examples/HydroTests/SodShock_BCC_3D/makeIC.py b/examples/HydroTests/SodShock_BCC_3D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..69f1bc506680d3f2f149c0fd7b75b069f9b00b64
--- /dev/null
+++ b/examples/HydroTests/SodShock_BCC_3D/makeIC.py
@@ -0,0 +1,119 @@
+###############################################################################
+ # 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/>.
+ #
+ ##############################################################################
+
+import h5py
+from numpy import *
+
+# Generates a swift IC file for the 3D Sod Shock in a periodic box
+
+# Parameters
+gamma = 5./3. # Gas adiabatic index
+x_min = -1.
+x_max = 1.
+rho_L = 1. # Density left state
+rho_R = 0.125 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
+fileName = "sodShock.hdf5"
+
+
+#---------------------------------------------------
+boxSize = (x_max - x_min)
+
+glass_L = h5py.File("glassCube_64.hdf5", "r")
+glass_R = h5py.File("glassCube_32.hdf5", "r")
+
+pos_L = glass_L["/PartType0/Coordinates"][:,:] * 0.5
+pos_R = glass_R["/PartType0/Coordinates"][:,:] * 0.5
+h_L = glass_L["/PartType0/SmoothingLength"][:] * 0.5
+h_R = glass_R["/PartType0/SmoothingLength"][:] * 0.5
+
+# Merge things
+aa = pos_L - array([0.5, 0., 0.])
+pos_LL = append(pos_L, pos_L + array([0.5, 0., 0.]), axis=0)
+pos_RR = append(pos_R, pos_R + array([0.5, 0., 0.]), axis=0)
+pos = append(pos_LL - array([1.0, 0., 0.]), pos_RR, axis=0)
+h_LL = append(h_L, h_L)
+h_RR = append(h_R, h_R)
+h = append(h_LL, h_RR)
+
+numPart_L = size(h_LL)
+numPart_R = size(h_RR)
+numPart = size(h)
+
+vol_L = 0.25
+vol_R = 0.25
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+u = zeros(numPart)
+
+for i in range(numPart):
+ x = pos[i,0]
+
+ if x < 0: #left
+ u[i] = P_L / (rho_L * (gamma - 1.))
+ m[i] = rho_L * vol_L / numPart_L
+ v[i,0] = v_L
+ else: #right
+ u[i] = P_R / (rho_R * (gamma - 1.))
+ m[i] = rho_R * vol_R / numPart_R
+ v[i,0] = v_R
+
+# Shift particles
+pos[:,0] -= x_min
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [boxSize, 0.5, 0.5]
+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
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 3
+
+#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.
+
+#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')
+
+
+file.close()
diff --git a/examples/HydroTests/SodShock_BCC_3D/plotSolution.py b/examples/HydroTests/SodShock_BCC_3D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..dc1b15c3eac862365d4422f95a14ffe1713f91a6
--- /dev/null
+++ b/examples/HydroTests/SodShock_BCC_3D/plotSolution.py
@@ -0,0 +1,221 @@
+###############################################################################
+# This file is part of the ANARCHY paper.
+# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+# 2019 Josh Borrow (joshua.boorrow@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/>.
+#
+################################################################################
+
+# Compares the swift result for the 2D spherical Sod shock with a high
+# resolution 2D reference result
+
+import sys
+import matplotlib
+
+matplotlib.use("Agg")
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy import stats
+
+from swiftsimio import load
+from analyticSolution import analytic
+
+snap = int(sys.argv[1])
+
+sim = load(f"sodshock_{snap:04d}.hdf5")
+
+# Set up plotting stuff
+try:
+ plt.style.use("mnras_durham")
+except:
+ rcParams = {
+ "font.serif": ["STIX", "Times New Roman", "Times"],
+ "font.family": ["serif"],
+ "mathtext.fontset": "stix",
+ "font.size": 8,
+ }
+ plt.rcParams.update(rcParams)
+
+
+# See analyticSolution for params.
+
+
+def get_data_dump(metadata):
+ """
+ Gets a big data dump from the SWIFT metadata
+ """
+
+ try:
+ viscosity = metadata.viscosity_info
+ except:
+ viscosity = "No info"
+
+ try:
+ diffusion = metadata.diffusion_info
+ except:
+ diffusion = "No info"
+
+ output = (
+ "$\\bf{SWIFT}$\n"
+ + metadata.code_info
+ + "\n\n"
+ + "$\\bf{Compiler}$\n"
+ + metadata.compiler_info
+ + "\n\n"
+ + "$\\bf{Hydrodynamics}$\n"
+ + metadata.hydro_info
+ + "\n\n"
+ + "$\\bf{Viscosity}$\n"
+ + viscosity
+ + "\n\n"
+ + "$\\bf{Diffusion}$\n"
+ + diffusion
+ )
+
+ return output
+
+
+# Read the simulation data
+time = sim.metadata.t.value
+
+data = dict(
+ x=sim.gas.coordinates.value[:, 0],
+ v=sim.gas.velocities.value[:, 0],
+ u=sim.gas.internal_energy.value,
+ S=sim.gas.entropy.value,
+ P=sim.gas.pressure.value,
+ rho=sim.gas.density.value,
+ y=sim.gas.coordinates.value[:, 1],
+ z=sim.gas.coordinates.value[:, 2],
+)
+
+# Try to add on the viscosity and diffusion.
+try:
+ data["visc"] = sim.gas.viscosity.value
+except:
+ pass
+
+try:
+ data["diff"] = 100.0 * sim.gas.diffusion.value
+except:
+ pass
+
+# Read in the "solution" data and calculate those that don't exist.
+
+ref = analytic(time=time)
+
+# We only want to plot this for the region that we actually have data for, hence the masking.
+mask = np.logical_and(ref["x"] < np.max(data["x"]), ref["x"] > np.min(data["x"]))
+ref = {k: v[mask] for k, v in ref.items()}
+
+# Now we can do the plotting.
+fig, ax = plt.subplots(2, 3, figsize=(6.974, 6.974 * (2.0 / 3.0)))
+ax = ax.flatten()
+
+# These are stored in priority order
+plot = dict(
+ v="Velocity ($v_x$)",
+ u="Internal Energy ($u$)",
+ rho=r"Density ($\rho$)",
+ P="Pressure ($P$)",
+ diff=r"100$\times$ Diffusion Coefficient ($\alpha_D$)",
+ visc=r"Viscosity Coefficient ($\alpha_V$)",
+ S="Entropy ($A$)",
+)
+
+log = dict(v=False, u=False, S=False, P=False, rho=False, visc=False, diff=False)
+ylim = dict(v=(-0.05, 1.0), diff=(0.0, None), visc=(0.0, None))
+
+current_axis = 0
+
+for key, label in plot.items():
+ if current_axis > 4:
+ break
+ else:
+ axis = ax[current_axis]
+
+ try:
+ if log[key]:
+ axis.semilogy()
+
+ # Raw data
+ axis.plot(
+ data["x"],
+ data[key],
+ ".",
+ color="C1",
+ markersize=0.5,
+ alpha=0.5,
+ rasterized=True,
+ markeredgecolor="none",
+ zorder=-1,
+ )
+
+ mask_noraster = np.logical_and.reduce([
+ data["y"] < 0.52,
+ data["y"] > 0.48,
+ data["z"] < 0.52,
+ data["z"] > 0.48
+ ])
+
+ axis.plot(
+ data["x"][mask_noraster],
+ data[key][mask_noraster],
+ ".",
+ color="C3",
+ rasterized=False,
+ markeredgecolor="none",
+ markersize=3,
+ zorder=0,
+ )
+
+ # Exact solution
+ try:
+ axis.plot(ref["x"], ref[key], c="C0", ls="dashed", zorder=1, lw=1)
+ except KeyError:
+ # No solution :(
+ pass
+
+ axis.set_xlabel("Position ($x$)", labelpad=0)
+ axis.set_ylabel(label, labelpad=0)
+
+ axis.set_xlim(0.6, 1.5)
+
+ try:
+ axis.set_ylim(*ylim[key])
+ except KeyError:
+ # No worries pal
+ pass
+
+ current_axis += 1
+ except KeyError:
+ # Mustn't have that data!
+ continue
+
+
+info_axis = ax[-1]
+
+info = get_data_dump(sim.metadata)
+
+info_axis.text(
+ 0.5, 0.45, info, ha="center", va="center", fontsize=5, transform=info_axis.transAxes
+)
+
+info_axis.axis("off")
+
+
+fig.tight_layout(pad=0.5)
+fig.savefig("SodShock.pdf", dpi=300)
diff --git a/examples/HydroTests/SodShock_BCC_3D/run.sh b/examples/HydroTests/SodShock_BCC_3D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..70cb3a97b1adf3944b9ed8965afabecfe280f7af
--- /dev/null
+++ b/examples/HydroTests/SodShock_BCC_3D/run.sh
@@ -0,0 +1,15 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e sodShock.hdf5 ]
+then
+ echo "Generating initial conditions for the Sod shock example..."
+ python makeGlass.py -n 64 -o glassCube_64.hdf5
+ python makeGlass.py -n 32 -o glassCube_32.hdf5
+ python makeIC.py
+fi
+
+# Run SWIFT
+../../swift --hydro --threads=4 sodShock.yml 2>&1 | tee output.log
+
+python plotSolution.py 1
diff --git a/examples/HydroTests/SodShock_BCC_3D/sodShock.yml b/examples/HydroTests/SodShock_BCC_3D/sodShock.yml
new file mode 100644
index 0000000000000000000000000000000000000000..373a70bf6b6782d7bdb4ed91417a13270f6521e6
--- /dev/null
+++ b/examples/HydroTests/SodShock_BCC_3D/sodShock.yml
@@ -0,0 +1,85 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1 # Grams
+ UnitLength_in_cgs: 1 # Centimeters
+ UnitVelocity_in_cgs: 1 # 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: 0.2 # The end time of the simulation (in internal units).
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sodShock # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.2 # Time difference between consecutive outputs (in internal units)
+ compression: 1
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
+# 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.
+ viscosity_alpha: 1.0
+ viscosity_alpha_max: 1.0
+ viscosity_alpha_min: 1.0
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./sodShock.hdf5 # The file to read
+ periodic: 1
+
+EAGLEChemistry: # Solar abundances
+ init_abundance_metal: 0.014
+ init_abundance_Hydrogen: 0.70649785
+ init_abundance_Helium: 0.28055534
+ init_abundance_Carbon: 2.0665436e-3
+ init_abundance_Nitrogen: 8.3562563e-4
+ init_abundance_Oxygen: 5.4926244e-3
+ init_abundance_Neon: 1.4144605e-3
+ init_abundance_Magnesium: 5.907064e-4
+ init_abundance_Silicon: 6.825874e-4
+ init_abundance_Iron: 1.1032152e-3
+
+EAGLECooling:
+ dir_name: ./coolingtables/
+ H_reion_z: 11.5
+ H_reion_eV_p_H: 2.0
+ He_reion_z_centre: 3.5
+ He_reion_z_sigma: 0.5
+ He_reion_eV_p_H: 2.0
+
+# EAGLE star formation parameters
+EAGLEStarFormation:
+ EOS_density_norm_H_p_cm3: 0.1 # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
+ EOS_temperature_norm_K: 8000 # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
+ EOS_gamma_effective: 1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
+ KS_normalisation: 1.515e-4 # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+ KS_exponent: 1.4 # The exponent of the Kennicutt-Schmidt law.
+ KS_min_over_density: 57.7 # The over-density above which star-formation is allowed.
+ KS_high_density_threshold_H_p_cm3: 1e3 # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+ KS_high_density_exponent: 2.0 # Slope of the Kennicut-Schmidt law above the high-density threshold.
+ KS_temperature_margin_dex: 0.5 # Logarithm base 10 of the maximal temperature difference above the EOS allowed to form stars.
+ threshold_norm_H_p_cm3: 0.1 # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+ threshold_Z0: 0.002 # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+ threshold_slope: -0.64 # Slope of the metal-dependant star formation threshold
+ threshold_max_density_H_p_cm3: 10.0 # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+
+# Parameters for the EAGLE "equation of state"
+EAGLEEntropyFloor:
+ Jeans_density_threshold_H_p_cm3: 0.1 # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Jeans_over_density_threshold: 10. # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
+ Jeans_temperature_norm_K: 8000 # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
+ Jeans_gamma_effective: 1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
+ Cool_density_threshold_H_p_cm3: 1e-5 # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+ Cool_over_density_threshold: 10. # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
+ Cool_temperature_norm_K: 8000 # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
+ Cool_gamma_effective: 1. # Slope the of the EAGLE Cool limiter entropy floor
diff --git a/examples/HydroTests/SquareTest_2D/README b/examples/HydroTests/SquareTest_2D/README
new file mode 100644
index 0000000000000000000000000000000000000000..fc07b9d9d17fc746752bd6eb47f03899b4b335a0
--- /dev/null
+++ b/examples/HydroTests/SquareTest_2D/README
@@ -0,0 +1,20 @@
+Square Test
+-----------
+
+This is a very challenging test that aims to figure out
+if contact discontinuities are properly handled. If there
+is residual surface tension, then the square will quickly
+become a sphere. Otherwise, it will remain a square. For
+more information see Hopkins' 2013 and 2015 papers.
+
+There are two initial condition generation files present.
+For the SWIFT method of finding an un-mass weighted number
+of particles in the kernel, it makes more sense to have
+different mass particles (makeICDifferentMasses.py). For
+comparison to previous methods, we also provide a script
+that creates initial conditions with a different density
+of particles, all with equal masses, in the square and
+outside of the square.
+
+If you do not have the swiftsimio library, you can use
+the plotSolutionLegacy.py to plot the solution.
diff --git a/examples/HydroTests/SquareTest_2D/makeIC.py b/examples/HydroTests/SquareTest_2D/makeIC.py
index 12a394873edf42f7ecfdf07c9795b62e3ad89745..bbb6837c019dde14d833588924901cb49e301b25 100644
--- a/examples/HydroTests/SquareTest_2D/makeIC.py
+++ b/examples/HydroTests/SquareTest_2D/makeIC.py
@@ -29,15 +29,15 @@ rho0 = 4 # Gas central density
rho1 = 1 # Gas outskirt density
P0 = 2.5 # Gas central pressure
P1 = 2.5 # Gas central pressure
-vx = 142.3 # Random velocity for all particles
-vy = -31.4
+vx = 0.0 # Random velocity for all particles
+vy = 0.0
fileOutputName = "square.hdf5"
#---------------------------------------------------
vol = 1.
numPart_out = L * L
-numPart_in = L * L * rho0 / rho1 / 4
+numPart_in = int(L * L * rho0 / rho1 / 4)
L_out = int(sqrt(numPart_out))
L_in = int(sqrt(numPart_in))
diff --git a/examples/HydroTests/SquareTest_2D/makeICDifferentMasses.py b/examples/HydroTests/SquareTest_2D/makeICDifferentMasses.py
new file mode 100644
index 0000000000000000000000000000000000000000..6c7af6d33ea020be773805d223095c99fbf78951
--- /dev/null
+++ b/examples/HydroTests/SquareTest_2D/makeICDifferentMasses.py
@@ -0,0 +1,110 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2019 Josh Borrow (joshua.borrow@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
+from numpy import *
+
+# Generates a swift IC file for the Square test in a periodic box
+
+# Parameters
+L = 2 * 64 # Number of particles on the side
+gamma = 5.0 / 3.0 # Gas adiabatic index
+rho0 = 4 # Gas central density
+rho1 = 1 # Gas outskirt density
+P0 = 2.5 # Gas central pressure
+P1 = 2.5 # Gas central pressure
+vx = 0.0 # Random velocity for all particles
+vy = 0.0
+fileOutputName = "square.hdf5"
+# ---------------------------------------------------
+
+vol = 1.0
+
+numPart = L * L
+
+pos_x = arange(0, 1, 1.0 / L)
+xv, yv = meshgrid(pos_x, pos_x)
+pos = zeros((numPart, 3), dtype=float)
+pos[:, 0] = xv.flatten()
+pos[:, 1] = yv.flatten()
+
+# Now we can get 2d masks!
+inside = logical_and.reduce([xv < 0.75, xv > 0.25, yv < 0.75, yv > 0.25])
+
+mass_in = rho0 / numPart
+mass_out = rho1 / numPart
+
+m = ones_like(xv) * mass_out
+m[inside] = mass_in
+m = m.flatten()
+
+h = ones_like(m) / L
+
+u_in = P0 / ((gamma - 1) * rho0)
+u_out = P1 / ((gamma - 1) * rho1)
+u = ones_like(xv) * u_out
+u[inside] = u_in
+u = u.flatten()
+
+vel = zeros_like(pos)
+
+ids = arange(numPart)
+
+
+# File
+fileOutput = h5py.File(fileOutputName, "w")
+
+# Header
+grp = fileOutput.create_group("/Header")
+grp.attrs["BoxSize"] = [vol, vol, 0.2]
+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
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["Dimension"] = 2
+
+# Units
+grp = fileOutput.create_group("/Units")
+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[()] = pos
+ds = grp.create_dataset("Velocities", (numPart, 3), "f")
+ds[()] = vel
+ds = grp.create_dataset("Masses", (numPart, 1), "f")
+ds[()] = m.reshape((numPart, 1))
+ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+ds[()] = h.reshape((numPart, 1))
+ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+ds[()] = u.reshape((numPart, 1))
+ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+ds[()] = ids.reshape((numPart, 1))
+
+fileOutput.close()
+
diff --git a/examples/HydroTests/SquareTest_2D/makeMovie.py b/examples/HydroTests/SquareTest_2D/makeMovie.py
new file mode 100644
index 0000000000000000000000000000000000000000..ccddca6e05b1103b8ece425c097f3ddec7bf2bb9
--- /dev/null
+++ b/examples/HydroTests/SquareTest_2D/makeMovie.py
@@ -0,0 +1,136 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2019 Josh Borrow (joshua.borrow@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/>.
+#
+##############################################################################
+
+from swiftsimio import load
+from swiftsimio.visualisation import project_gas_pixel_grid
+from matplotlib.animation import FuncAnimation
+import matplotlib.pyplot as plt
+
+from numpy import max, min
+
+try:
+ # Try and load this, otherwise we're stuck with serial
+ from p_tqdm import p_map
+
+ map = p_map
+except:
+ print("Try installing the p_tqdm package to make movie frames in parallel")
+ pass
+
+
+def load_and_make_image(number, res, property):
+ filename = "square_{:04d}.hdf5".format(number)
+ data = load(filename)
+ image = project_gas_pixel_grid(data, res, property)
+ image_none = project_gas_pixel_grid(data, res, None)
+
+ return image / image_none
+
+
+def create_movie(filename, start, stop, resolution, property, output_filename):
+ """
+ Creates the movie with:
+
+ snapshots named {filename}_{start}.hdf5 to {filename}_{stop}.hdf5
+ at {resolution} x {resolution} pixel size and smoothing the given
+ {property} and outputting to {output_filename}.
+ """
+
+ def baked_in_load(n):
+ return load_and_make_image(n, resolution, property)
+
+ # Make frames in parallel (reading also parallel!)
+ frames = map(baked_in_load, list(range(start, stop)))
+
+ vmax = max(list(map(max, frames)))
+ vmin = min(list(map(min, frames)))
+
+ fig, ax = plt.subplots(figsize=(1, 1))
+ ax.axis("off")
+ fig.subplots_adjust(0, 0, 1, 1)
+
+ image = ax.imshow(frames[0], origin="lower", vmax=vmax, vmin=vmin)
+
+ def frame(n):
+ image.set_array(frames[n])
+
+ return (image,)
+
+ animation = FuncAnimation(fig, frame, range(0, 40), interval=40)
+ animation.save(output_filename, dpi=resolution)
+
+
+if __name__ == "__main__":
+ from argparse import ArgumentParser
+
+ parser = ArgumentParser(description="""Creates a movie of the whole box.""")
+
+ parser.add_argument(
+ "-s",
+ "--snapshot",
+ help="Snapshot name. Default: square",
+ type=str,
+ default="square",
+ )
+
+ parser.add_argument(
+ "-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
+ )
+
+ parser.add_argument(
+ "-o",
+ "--output",
+ help="Output filename. Default: EnergyMovie.mp4",
+ type=str,
+ default="EnergyMovie.mp4",
+ )
+
+ parser.add_argument(
+ "-p",
+ "--property",
+ help="(swiftsimio) Property to plot. Default: internal_energy",
+ type=str,
+ default="internal_energy",
+ )
+
+ parser.add_argument(
+ "-r",
+ "--resolution",
+ help="Resolution to make the movie at. Default: 512",
+ type=int,
+ default=512,
+ )
+
+ vars = parser.parse_args()
+
+ create_movie(
+ filename=vars.snapshot,
+ start=vars.initial,
+ stop=vars.final,
+ resolution=vars.resolution,
+ property=vars.property,
+ output_filename=vars.output,
+ )
+
+ exit(0)
+
diff --git a/examples/HydroTests/SquareTest_2D/plotSolution.py b/examples/HydroTests/SquareTest_2D/plotSolution.py
index f182b4d7437348d29065b51df79e5334aa26f9a4..ae351abb47eb5bf63f71b55b33a060067c8a6a01 100644
--- a/examples/HydroTests/SquareTest_2D/plotSolution.py
+++ b/examples/HydroTests/SquareTest_2D/plotSolution.py
@@ -1,186 +1,171 @@
###############################################################################
- # 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/>.
- #
- ##############################################################################
-
-# Computes the analytical solution of the square test
-
-# Parameters
-gas_gamma = 5./3. # Gas adiabatic index
-gamma = 5./3. # Gas adiabatic index
-rho0 = 4 # Gas central density
-rho1 = 1 # Gas outskirt density
-P0 = 2.5 # Gas central pressure
-P1 = 2.5 # Gas central pressure
-vx = 142.3 # Random velocity for all particles
-vy = -31.4
-
-# ---------------------------------------------------------------
-# Don't touch anything after this.
-# ---------------------------------------------------------------
-
+# This file is part of SWIFT.
+# Copyright (c) 2019 Josh Borrow (joshua.borrow@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/>.
+#
+##############################################################################
+"""
+Plots the solution of the square test in a smoothed way using SWIFTsimIO's
+smoothed plotting.
+"""
+
+import sys
import matplotlib
+
matplotlib.use("Agg")
-from pylab import *
-import h5py
-
-# 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']})
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy import stats
+
+from swiftsimio import load
+from swiftsimio.visualisation import project_gas_pixel_grid
snap = int(sys.argv[1])
-# Read the simulation data
-sim = h5py.File("square_%04d.hdf5"%snap, "r")
-boxSize = sim["/Header"].attrs["BoxSize"][0]
-time = sim["/Header"].attrs["Time"][0]
-scheme = sim["/HydroScheme"].attrs["Scheme"]
-kernel = sim["/HydroScheme"].attrs["Kernel function"]
-neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
-eta = sim["/HydroScheme"].attrs["Kernel eta"]
-git = sim["Code"].attrs["Git Revision"]
-
-# Analytical soltion
-centre_x = 0.5 + time * vx
-centre_y = 0.5 + time * vy
-
-while centre_x > 1.:
- centre_x -= 1.
-while centre_x < 0.:
- centre_x += 1.
-while centre_y > 1.:
- centre_y -= 1.
-while centre_y < 0.:
- centre_y += 1.
-
-pos = sim["/PartType0/Coordinates"][:,:]
-vel = sim["/PartType0/Velocities"][:,:]
-v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
-rho = sim["/PartType0/Density"][:]
-u = sim["/PartType0/InternalEnergy"][:]
-S = sim["/PartType0/Entropy"][:]
-P = sim["/PartType0/Pressure"][:]
-x = pos[:,0] - centre_x
-y = pos[:,1] - centre_y
-
-# Box wrapping
-x[x>0.5] -= 1.
-x[x<-0.5] += 1.
-y[y>0.5] -= 1.
-y[y<-0.5] += 1.
-
-# Azimuthal velocity profile -----------------------------
-subplot(231)
-scatter(x, y, c=v_norm, cmap="PuBu", edgecolors='face', s=4, vmin=-1., vmax=1.)
-text(0.47, 0.47, "${\\rm{Velocity~norm}}$", ha="right", va="top", backgroundcolor="w")
-plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
-plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=-7)
-xlim(-0.5, 0.5)
-ylim(-0.5, 0.5)
-
-# Radial density profile --------------------------------
-subplot(232)
-scatter(x, y, c=rho, cmap="PuBu", edgecolors='face', s=4, vmin=0., vmax=4.)
-text(0.47, 0.47, "${\\rm{Density}}$", ha="right", va="top", backgroundcolor="w")
-plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
-plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=-7)
-xlim(-0.5, 0.5)
-ylim(-0.5, 0.5)
-
-# Radial pressure profile --------------------------------
-subplot(233)
-scatter(x, y, c=P, cmap="PuBu", edgecolors='face', s=4, vmin=2, vmax=4)
-text(0.47, 0.47, "${\\rm{Pressure}}$", ha="right", va="top", backgroundcolor="w")
-plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
-plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=-7)
-xlim(-0.5, 0.5)
-ylim(-0.5, 0.5)
-
-# Internal energy profile --------------------------------
-subplot(234)
-scatter(x, y, c=u, cmap="PuBu", edgecolors='face', s=4, vmin=0.5, vmax=4.)
-text(0.47, 0.47, "${\\rm{Internal~energy}}$", ha="right", va="top", backgroundcolor="w")
-plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
-plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=-7)
-xlim(-0.5, 0.5)
-ylim(-0.5, 0.5)
-
-# Radial entropy profile --------------------------------
-subplot(235)
-scatter(x, y, c=S, cmap="PuBu", edgecolors='face', s=4, vmin=0., vmax=3.)
-text(0.47, 0.47, "${\\rm{Entropy}}$", ha="right", va="top", backgroundcolor="w")
-plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
-plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
-xlabel("${\\rm{Position}}~x$", labelpad=0)
-ylabel("${\\rm{Position}}~y$", labelpad=-7)
-xlim(-0.5, 0.5)
-ylim(-0.5, 0.5)
-
-
-# Information -------------------------------------
-subplot(236, frameon=False)
-
-text(-0.49, 0.9, "Square test with $\\gamma=%.3f$ at $t=%.2f$"%(gas_gamma,time), fontsize=10)
-text(-0.49, 0.8, "Centre:~~~ $(P, \\rho) = (%.3f, %.3f)$"%(P0, rho0), fontsize=10)
-text(-0.49, 0.7, "Outskirts: $(P, \\rho) = (%.3f, %.3f)$"%(P1, rho1), 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("SquareTest.png", dpi=200)
+sim = load(f"square_{snap:04d}.hdf5")
+resolution = 512
+
+# First create a grid that gets the particle density so we can divide it out later
+unweighted_grid = project_gas_pixel_grid(sim, 512, None)
+
+# Set up plotting stuff
+try:
+ plt.style.use("mnras_durham")
+except:
+ rcParams = {
+ "font.serif": ["STIX", "Times New Roman", "Times"],
+ "font.family": ["serif"],
+ "mathtext.fontset": "stix",
+ "font.size": 8,
+ }
+ plt.rcParams.update(rcParams)
+
+
+def get_data_dump(metadata):
+ """
+ Gets a big data dump from the SWIFT metadata
+ """
+
+ try:
+ viscosity = metadata.viscosity_info
+ except:
+ viscosity = "No info"
+
+ try:
+ diffusion = metadata.diffusion_info
+ except:
+ diffusion = "No info"
+
+ output = (
+ "$\\bf{SWIFT}$\n"
+ + metadata.code_info
+ + "\n\n"
+ + "$\\bf{Compiler}$\n"
+ + metadata.compiler_info
+ + "\n\n"
+ + "$\\bf{Hydrodynamics}$\n"
+ + metadata.hydro_info
+ + "\n\n"
+ + "$\\bf{Viscosity}$\n"
+ + viscosity
+ + "\n\n"
+ + "$\\bf{Diffusion}$\n"
+ + diffusion
+ )
+
+ return output
+
+
+# Now we can do the plotting.
+fig, ax = plt.subplots(2, 3, figsize=(6.974, 6.974 * (2.0 / 3.0)))
+ax = ax.flatten()
+
+# These are stored in priority order
+plot = dict(
+ internal_energy="Internal Energy ($u$)",
+ density=r"Density ($\rho$)",
+ pressure="Pressure ($P$)",
+ entropy="Entropy ($A$)",
+ diffusion=r"Diffusion ($\alpha$)",
+)
+
+current_axis = 0
+
+for key, label in plot.items():
+ if current_axis > 4:
+ break
+ else:
+ axis = ax[current_axis]
+
+ try:
+ # Raw data
+ try:
+ grid = (
+ project_gas_pixel_grid(sim, resolution=resolution, project=key)
+ / unweighted_grid
+ )
+ axis.imshow(grid, origin="lower", extent=[0, 1, 0, 1], cmap="cividis")
+ except:
+ continue
+
+ # Exact solution, a square!
+ axis.plot(
+ [0.25, 0.75, 0.75, 0.25, 0.25],
+ [0.25, 0.25, 0.75, 0.75, 0.25],
+ linestyle="dashed",
+ color="white",
+ )
+
+ circle = plt.Circle(
+ (0.8, 0.8),
+ radius=np.sqrt(1.0 / sim.metadata.n_gas) * 1.25,
+ linestyle="solid",
+ color="white",
+ fill=False,
+ )
+
+ axis.add_artist(circle)
+
+ axis.tick_params(
+ axis="both",
+ which="both",
+ labelbottom=False,
+ labelleft=False,
+ bottom=False,
+ left=False,
+ )
+
+ axis.set_title(label)
+
+ current_axis += 1
+ except KeyError:
+ # Mustn't have that data!
+ continue
+
+
+info_axis = ax[-1]
+
+info = get_data_dump(sim.metadata)
+
+info_axis.text(
+ 0.5, 0.45, info, ha="center", va="center", fontsize=5, transform=info_axis.transAxes
+)
+
+info_axis.axis("off")
+
+
+fig.tight_layout(pad=1.0)
+fig.savefig("SquareTest.pdf", dpi=300)
diff --git a/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py b/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py
new file mode 100644
index 0000000000000000000000000000000000000000..956da800c9096232d2e82cf4cff4c780672e0a8f
--- /dev/null
+++ b/examples/HydroTests/SquareTest_2D/plotSolutionLegacy.py
@@ -0,0 +1,186 @@
+###############################################################################
+ # 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/>.
+ #
+ ##############################################################################
+
+# Computes the analytical solution of the square test
+
+# Parameters
+gas_gamma = 5./3. # Gas adiabatic index
+gamma = 5./3. # Gas adiabatic index
+rho0 = 4 # Gas central density
+rho1 = 1 # Gas outskirt density
+P0 = 2.5 # Gas central pressure
+P1 = 2.5 # Gas central pressure
+vx = 0.0 # Random velocity for all particles
+vy = 0.0
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# 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']})
+
+snap = int(sys.argv[1])
+
+# Read the simulation data
+sim = h5py.File("square_%04d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+
+# Analytical soltion
+centre_x = 0.5 + time * vx
+centre_y = 0.5 + time * vy
+
+while centre_x > 1.:
+ centre_x -= 1.
+while centre_x < 0.:
+ centre_x += 1.
+while centre_y > 1.:
+ centre_y -= 1.
+while centre_y < 0.:
+ centre_y += 1.
+
+pos = sim["/PartType0/Coordinates"][:,:]
+vel = sim["/PartType0/Velocities"][:,:]
+v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
+rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+x = pos[:,0] - centre_x
+y = pos[:,1] - centre_y
+
+# Box wrapping
+x[x>0.5] -= 1.
+x[x<-0.5] += 1.
+y[y>0.5] -= 1.
+y[y<-0.5] += 1.
+
+# Azimuthal velocity profile -----------------------------
+subplot(231)
+scatter(x, y, c=v_norm, cmap="PuBu", edgecolors='face', s=4, vmin=-1., vmax=1.)
+text(0.47, 0.47, "${\\rm{Velocity~norm}}$", ha="right", va="top", backgroundcolor="w")
+plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
+plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=-7)
+xlim(-0.5, 0.5)
+ylim(-0.5, 0.5)
+
+# Radial density profile --------------------------------
+subplot(232)
+scatter(x, y, c=rho, cmap="PuBu", edgecolors='face', s=4, vmin=0., vmax=4.)
+text(0.47, 0.47, "${\\rm{Density}}$", ha="right", va="top", backgroundcolor="w")
+plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
+plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=-7)
+xlim(-0.5, 0.5)
+ylim(-0.5, 0.5)
+
+# Radial pressure profile --------------------------------
+subplot(233)
+scatter(x, y, c=P, cmap="PuBu", edgecolors='face', s=4, vmin=2, vmax=4)
+text(0.47, 0.47, "${\\rm{Pressure}}$", ha="right", va="top", backgroundcolor="w")
+plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
+plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=-7)
+xlim(-0.5, 0.5)
+ylim(-0.5, 0.5)
+
+# Internal energy profile --------------------------------
+subplot(234)
+scatter(x, y, c=u, cmap="PuBu", edgecolors='face', s=4, vmin=0.5, vmax=4.)
+text(0.47, 0.47, "${\\rm{Internal~energy}}$", ha="right", va="top", backgroundcolor="w")
+plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
+plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=-7)
+xlim(-0.5, 0.5)
+ylim(-0.5, 0.5)
+
+# Radial entropy profile --------------------------------
+subplot(235)
+scatter(x, y, c=S, cmap="PuBu", edgecolors='face', s=4, vmin=0., vmax=3.)
+text(0.47, 0.47, "${\\rm{Entropy}}$", ha="right", va="top", backgroundcolor="w")
+plot([-0.25, 0.25], [0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, 0.25], [-0.25, -0.25], '--', color='k', alpha=0.8, lw=2)
+plot([0.25, 0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+plot([-0.25, -0.25], [-0.25, 0.25], '--', color='k', alpha=0.8, lw=2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=-7)
+xlim(-0.5, 0.5)
+ylim(-0.5, 0.5)
+
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Square test with $\\gamma=%.3f$ at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "Centre:~~~ $(P, \\rho) = (%.3f, %.3f)$"%(P0, rho0), fontsize=10)
+text(-0.49, 0.7, "Outskirts: $(P, \\rho) = (%.3f, %.3f)$"%(P1, rho1), 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("SquareTest.png", dpi=200)
diff --git a/examples/HydroTests/SquareTest_2D/run.sh b/examples/HydroTests/SquareTest_2D/run.sh
index dae0d742e706d73f3a5efc26a9d82ac59c883757..f5e0160627f7417af4f7bf38e50c03214e2f3c90 100755
--- a/examples/HydroTests/SquareTest_2D/run.sh
+++ b/examples/HydroTests/SquareTest_2D/run.sh
@@ -4,11 +4,12 @@
if [ ! -e square.hdf5 ]
then
echo "Generating initial conditions for the square test ..."
- python makeIC.py
+ python makeICDifferentMasses.py
fi
# Run SWIFT
-../../swift --hydro --threads=1 square.yml 2>&1 | tee output.log
+../../swift --hydro --threads=4 square.yml 2>&1 | tee output.log
# Plot the solution
-python plotSolution.py 5
+python plotSolution.py 40
+python makeMovie.py