Merge branch 'add-sodshock-bcc' into 'master'

Adds the SodShock BCC test and updates SquareTest

See merge request swift/swiftsim!840

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")

examples/HydroTests/SodShock_BCC_3D/README

+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.

examples/HydroTests/SodShock_BCC_3D/analyticSolution.py

+###############################################################################
+# 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)

examples/HydroTests/SodShock_BCC_3D/makeGlass.py

+"""
+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)
+

examples/HydroTests/SodShock_BCC_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/>.
+ # 
+ ##############################################################################
+
+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()

examples/HydroTests/SodShock_BCC_3D/plotSolution.py

+###############################################################################
+# 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)

examples/HydroTests/SodShock_BCC_3D/run.sh

+#!/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

examples/HydroTests/SodShock_BCC_3D/sodShock.yml

+# 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

examples/HydroTests/SquareTest_2D/README

+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.

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))

examples/HydroTests/SquareTest_2D/makeICDifferentMasses.py

+###############################################################################
+# 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()
+

examples/HydroTests/SquareTest_2D/makeMovie.py

+###############################################################################
+# 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)
+

examples/HydroTests/SquareTest_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/>.
- # 
- ##############################################################################
-
-# 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)

examples/HydroTests/SquareTest_2D/plotSolutionLegacy.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/>.
+ # 
+ ##############################################################################
+
+# 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)

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