Added reference result for interacting blast waves. Added 2D and 3D spherical...

Added reference result for interacting blast waves. Added 2D and 3D spherical vacuum expansion test. Disabled pathological particle configuration errors in GIZMO to make 3D vacuum test work.

examples/EvrardCollapse_3D/plotSolution.py

@@ -110,7 +110,7 @@ semilogx(x, -v, '.', color='r', ms=0.2)
 semilogx(ref[:,0], ref[:,2], "k--", alpha=0.8, lw=1.2)
 errorbar(x_bin, -v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
+ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
 xlim(1.e-3, 2.)
 ylim(-1.7, 0.1)
 
@@ -157,7 +157,7 @@ ylim(0., 0.25)
 # Information -------------------------------------
 subplot(236, frameon=False)
 
-text(-0.49, 0.9, "Evrard collapse with  $\\gamma=%.3f$ in 3D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.9, "Evrard collapse with $\\gamma=%.3f$ in 3D\nat $t=%.2f$"%(gas_gamma,time), 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)
@@ -168,5 +168,5 @@ ylim(0, 1)
 xticks([])
 yticks([])
 
-
+tight_layout()
 savefig("EvrardCollapse.png", dpi=200)

examples/InteractingBlastWaves_1D/getReference.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/interactingBlastWaves1D_exact.txt

examples/InteractingBlastWaves_1D/plotSolution.py

@@ -17,21 +17,116 @@
 #
 ##############################################################################
 
-import h5py
+import numpy as np
 import matplotlib
 matplotlib.use("Agg")
 import pylab as pl
+import h5py
 import sys
 
+# Parameters
+gamma = 1.4 # Polytropic index
+
+# 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
+}
+pl.rcParams.update(params)
+pl.rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+# Read the snapshot index from the command line argument
 snap = int(sys.argv[1])
 
-sim = h5py.File("interactingBlastWaves_%04d.hdf5" % snap, "r")
-coords = sim["/PartType0/Coordinates"]
-rho = sim["/PartType0/Density"]
+# Open the file and read the relevant data
+file = h5py.File("interactingBlastWaves_{0:04d}.hdf5".format(snap), "r")
+x = file["/PartType0/Coordinates"][:,0]
+rho = file["/PartType0/Density"]
+v = file["/PartType0/Velocities"][:,0]
+u = file["/PartType0/InternalEnergy"]
+S = file["/PartType0/Entropy"]
+P = file["/PartType0/Pressure"]
+time = file["/Header"].attrs["Time"][0]
+
+scheme = file["/HydroScheme"].attrs["Scheme"]
+kernel = file["/HydroScheme"].attrs["Kernel function"]
+neighbours = file["/HydroScheme"].attrs["Kernel target N_ngb"][0]
+eta = file["/HydroScheme"].attrs["Kernel eta"][0]
+git = file["Code"].attrs["Git Revision"]
+
+ref = np.loadtxt("interactingBlastWaves1D_exact.txt")
+
+# Plot the interesting quantities
+fig, ax = pl.subplots(2, 3)
+
+# Velocity profile
+ax[0][0].plot(x, v, "r.", markersize = 4.)
+ax[0][0].plot(ref[:,0], ref[:,2], "k--", alpha = 0.8, linewidth = 1.2)
+ax[0][0].set_xlabel("${\\rm{Position}}~x$", labelpad = 0)
+ax[0][0].set_ylabel("${\\rm{Velocity}}~v_x$", labelpad = 0)
+ax[0][0].set_xlim(0., 1.)
+ax[0][0].set_ylim(-1., 15.)
+
+# Density profile
+ax[0][1].plot(x, rho, "r.", markersize = 4.)
+ax[0][1].plot(ref[:,0], ref[:,1], "k--", alpha = 0.8, linewidth = 1.2)
+ax[0][1].set_xlabel("${\\rm{Position}}~x$", labelpad = 0)
+ax[0][1].set_ylabel("${\\rm{Density}}~\\rho$", labelpad = 0)
+ax[0][1].set_xlim(0., 1.)
+
+# Pressure profile
+ax[0][2].plot(x, P, "r.", markersize = 4.)
+ax[0][2].plot(ref[:,0], ref[:,3], "k--", alpha = 0.8, linewidth = 1.2)
+ax[0][2].set_xlabel("${\\rm{Position}}~x$", labelpad = 0)
+ax[0][2].set_ylabel("${\\rm{Pressure}}~P$", labelpad = 0)
+ax[0][2].set_xlim(0., 1.)
+
+# Internal energy profile
+ax[1][0].plot(x, u, "r.", markersize = 4.)
+ax[1][0].plot(ref[:,0], ref[:,3] / ref[:,1] / (gamma - 1.), "k--", alpha = 0.8,
+              linewidth = 1.2)
+ax[1][0].set_xlabel("${\\rm{Position}}~x$", labelpad = 0)
+ax[1][0].set_ylabel("${\\rm{Internal~Energy}}~u$", labelpad = 0)
+ax[1][0].set_xlim(0., 1.)
+
+# Entropy profile
+ax[1][1].plot(x, S, "r.", markersize = 4.)
+ax[1][1].plot(ref[:,0], ref[:,3] / ref[:,1]**gamma, "k--", alpha = 0.8,
+              linewidth = 1.2)
+ax[1][1].set_xlabel("${\\rm{Position}}~x$", labelpad = 0)
+ax[1][1].set_ylabel("${\\rm{Entropy}}~S$", labelpad = 0)
+ax[1][1].set_xlim(0., 1.)
+
+# Run information
+ax[1][2].set_frame_on(False)
+ax[1][2].text(-0.49, 0.9,
+  "Interacting blast waves test\nwith $\\gamma={0:.3f}$ in 1D at $t = {1:.2f}$".format(
+    gamma, time), fontsize = 10)
+ax[1][2].plot([-0.49, 0.1], [0.62, 0.62], "k-", lw = 1)
+ax[1][2].text(-0.49, 0.5, "$\\textsc{{Swift}}$ {0}".format(git), fontsize = 10)
+ax[1][2].text(-0.49, 0.4, scheme, fontsize = 10)
+ax[1][2].text(-0.49, 0.3, kernel, fontsize = 10)
+ax[1][2].text(-0.49, 0.2,
+  "${0:.2f}$ neighbours ($\\eta={1:.3f}$)".format(neighbours, eta),
+  fontsize = 10)
+ax[1][2].set_xlim(-0.5, 0.5)
+ax[1][2].set_ylim(0., 1.)
+ax[1][2].set_xticks([])
+ax[1][2].set_yticks([])
 
-pl.xlabel("$x$")
-pl.ylabel("$\\rho{}$")
-pl.xlim(0.4, 1.)
-pl.ylim(0., 7.)
-pl.plot(coords[:,0], rho, "k.")
-pl.savefig("InteractingBlastWaves.png")
+pl.tight_layout()
+pl.savefig("InteractingBlastWaves.png", dpi = 200)

examples/InteractingBlastWaves_1D/run.sh

 #!/bin/bash
 
- # Generate the initial conditions if they are not present.
+# Generate the initial conditions if they are not present.
 if [ ! -e interactingBlastWaves.hdf5 ]
 then
     echo "Generating initial conditions for the Sedov blast example..."
@@ -10,5 +10,12 @@ fi
 # Run SWIFT
 ../swift -s -t 1 interactingBlastWaves.yml 2>&1 | tee output.log
 
+# Get the high resolution reference solution if not present.
+if [ ! -e interactingBlastWaves1D_exact.txt ]
+then
+    echo "Fetching reference solution for the interacting blast waves test..."
+    ./getReference.sh
+fi
+
 # Plot the solution
 python plotSolution.py 4

examples/SodShockSpherical_2D/plotSolution.py

@@ -109,7 +109,7 @@ plot(x, v, '.', color='r', ms=0.2)
 plot(ref[:,0], ref[:,2], "k--", alpha=0.8, lw=1.2)
 errorbar(x_bin, v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
+ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
 
 # Density profile --------------------------------
 subplot(232)

examples/SodShockSpherical_3D/plotSolution.py

@@ -110,7 +110,7 @@ plot(x, v, '.', color='r', ms=0.2)
 plot(ref[:,0], ref[:,2], "k--", alpha=0.8, lw=1.2)
 errorbar(x_bin, v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
-ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
+ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
 
 # Density profile --------------------------------
 subplot(232)

examples/VacuumSpherical_2D/getGlass.sh

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

examples/VacuumSpherical_2D/getReference.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/vacuumSpherical2D_exact.txt

examples/VacuumSpherical_2D/makeIC.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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 numpy as np
+import h5py
+
+# Generates an overdensity within a vacuum to test the vacuum resolving
+# capabilities of the code
+
+# Parameters
+gamma = 5. / 3.    # Gas adiabatic index
+
+fileName = "vacuum.hdf5" 
+
+#---------------------------------------------------
+glass = h5py.File("glassPlane_128.hdf5", "r")
+
+# Read particle positions and h from the glass
+pos = glass["/PartType0/Coordinates"][:,:]
+h = glass["/PartType0/SmoothingLength"][:] * 0.3
+
+# Make 4 copies of the glass to have more particles
+pos *= 0.5
+h *= 0.5
+pos = np.append(pos, pos + np.array([0.5, 0., 0.]), axis = 0)
+h = np.append(h, h)
+pos = np.append(pos, pos + np.array([0., 0.5, 0.]), axis = 0)
+h = np.append(h, h)
+
+radius = np.sqrt((pos[:,0] - 0.5)**2 + (pos[:,1] - 0.5)**2)
+index = radius < 0.25
+pos = pos[index]
+h = h[index]
+
+numPart = len(h)
+vol = np.pi * 0.25**2
+
+# Generate extra arrays
+v = np.zeros((numPart, 3))
+ids = np.linspace(1, numPart, numPart)
+m = np.zeros(numPart)
+u = np.zeros(numPart)
+
+m[:] = 1. * vol / numPart
+u[:] = 1. / (1. * (gamma - 1.))
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [1., 1., 1.]
+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"] = 2
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#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/VacuumSpherical_2D/plotSolution.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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 numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import pylab as pl
+import h5py
+import sys
+import scipy.stats as stats
+
+# Parameters
+gamma = 5. / 3. # Polytropic index
+rhoL = 1.       # Initial density in the non vacuum state
+vL = 0.         # Initial velocity in the non vacuum state
+PL = 1.         # Initial pressure in the non vacuum state
+rhoR = 0.       # Initial vacuum density
+vR = 0.         # Initial vacuum velocity
+PR = 0.         # Initial vacuum pressure
+
+# 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
+}
+pl.rcParams.update(params)
+pl.rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+# Read the snapshot index from the command line argument
+snap = int(sys.argv[1])
+
+# Open the file and read the relevant data
+file = h5py.File("vacuum_{0:04d}.hdf5".format(snap), "r")
+coords = file["/PartType0/Coordinates"]
+x = np.sqrt((coords[:,0] - 0.5)**2 + (coords[:,1] - 0.5)**2)
+rho = file["/PartType0/Density"][:]
+vels = file["/PartType0/Velocities"]
+v = np.sqrt(vels[:,0]**2 + vels[:,1]**2)
+u = file["/PartType0/InternalEnergy"][:]
+S = file["/PartType0/Entropy"][:]
+P = file["/PartType0/Pressure"][:]
+time = file["/Header"].attrs["Time"][0]
+
+scheme = file["/HydroScheme"].attrs["Scheme"]
+kernel = file["/HydroScheme"].attrs["Kernel function"]
+neighbours = file["/HydroScheme"].attrs["Kernel target N_ngb"][0]
+eta = file["/HydroScheme"].attrs["Kernel eta"][0]
+git = file["Code"].attrs["Git Revision"]
+
+# Bin the data values
+# We let scipy choose the bins and then reuse them for all other quantities
+rho_bin, x_bin_edge, _ = \
+  stats.binned_statistic(x, rho, statistic = "mean", bins = 50)
+rho2_bin, _, _ = \
+  stats.binned_statistic(x, rho**2, statistic = "mean", bins = x_bin_edge)
+rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
+
+v_bin, _, _ = \
+  stats.binned_statistic(x, v, statistic = "mean", bins = x_bin_edge)
+v2_bin, _, _ = \
+  stats.binned_statistic(x, v**2, statistic = "mean", bins = x_bin_edge)
+v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
+
+P_bin, _, _ = \
+  stats.binned_statistic(x, P, statistic = "mean", bins = x_bin_edge)
+P2_bin, _, _ = \
+  stats.binned_statistic(x, P**2, statistic = "mean", bins = x_bin_edge)
+P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
+
+u_bin, _, _ = \
+  stats.binned_statistic(x, u, statistic = "mean", bins = x_bin_edge)
+u2_bin, _, _ = \
+  stats.binned_statistic(x, u**2, statistic = "mean", bins = x_bin_edge)
+u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
+
+S_bin, _, _ = \
+  stats.binned_statistic(x, S, statistic = "mean", bins = x_bin_edge)
+S2_bin, _, _ = \
+  stats.binned_statistic(x, S**2, statistic = "mean", bins = x_bin_edge)
+S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
+
+x_bin = 0.5 * (x_bin_edge[1:] + x_bin_edge[:-1])
+
+ref = np.loadtxt("vacuumSpherical2D_exact.txt")
+
+# Plot the interesting quantities
+fig, ax = pl.subplots(2, 3)
+
+# Velocity profile
+ax[0][0].plot(x, v, "r.", markersize = 0.2)
+ax[0][0].plot(ref[:,0], ref[:,2], "k--", alpha = 0.8, linewidth = 1.2)
+ax[0][0].errorbar(x_bin, v_bin, yerr = v_sigma_bin, fmt = ".",
+  markersize = 8., color = "b", linewidth = 1.2)
+ax[0][0].set_xlabel("${\\rm{Radius}}~r$", labelpad = 0)
+ax[0][0].set_ylabel("${\\rm{Velocity}}~v_r$", labelpad = 0)
+ax[0][0].set_xlim(0., 0.4)
+ax[0][0].set_ylim(-0.1, 3.2)
+
+# Density profile
+ax[0][1].plot(x, rho, "r.", markersize = 0.2)
+ax[0][1].plot(ref[:,0], ref[:,1], "k--", alpha = 0.8, linewidth = 1.2)
+ax[0][1].errorbar(x_bin, rho_bin, yerr = rho_sigma_bin, fmt = ".",
+  markersize = 8., color = "b", linewidth = 1.2)
+ax[0][1].set_xlabel("${\\rm{Radius}}~r$", labelpad = 0)
+ax[0][1].set_ylabel("${\\rm{Density}}~\\rho$", labelpad = 0)
+ax[0][1].set_xlim(0., 0.4)
+
+# Pressure profile
+ax[0][2].plot(x, P, "r.", markersize = 0.2)
+ax[0][2].plot(ref[:,0], ref[:,3], "k--", alpha = 0.8, linewidth = 1.2)
+ax[0][2].errorbar(x_bin, P_bin, yerr = P_sigma_bin, fmt = ".",
+  markersize = 8., color = "b", linewidth = 1.2)
+ax[0][2].set_xlabel("${\\rm{Radius}}~r$", labelpad = 0)
+ax[0][2].set_ylabel("${\\rm{Pressure}}~P$", labelpad = 0)
+ax[0][2].set_xlim(0., 0.4)
+
+# Internal energy profile
+ax[1][0].plot(x, u, "r.", markersize = 0.2)
+ax[1][0].plot(ref[:,0], ref[:,3] / ref[:,1] / (gamma - 1.), "k--", alpha = 0.8,
+              linewidth = 1.2)
+ax[1][0].errorbar(x_bin, u_bin, yerr = u_sigma_bin, fmt = ".",
+  markersize = 8., color = "b", linewidth = 1.2)
+ax[1][0].set_xlabel("${\\rm{Radius}}~r$", labelpad = 0)
+ax[1][0].set_ylabel("${\\rm{Internal~Energy}}~u$", labelpad = 0)
+ax[1][0].set_xlim(0., 0.4)
+
+# Entropy profile
+ax[1][1].plot(x, S, "r.", markersize = 0.2)
+ax[1][1].plot(ref[:,0], ref[:,3] / ref[:,1]**gamma, "k--", alpha = 0.8,
+              linewidth = 1.2)
+ax[1][1].errorbar(x_bin, S_bin, yerr = S_sigma_bin, fmt = ".",
+  markersize = 8., color = "b", linewidth = 1.2)
+ax[1][1].set_xlabel("${\\rm{Radius}}~r$", labelpad = 0)
+ax[1][1].set_ylabel("${\\rm{Entropy}}~S$", labelpad = 0)
+ax[1][1].set_xlim(0., 0.4)
+ax[1][1].set_ylim(0., 4.)
+
+# Run information
+ax[1][2].set_frame_on(False)
+ax[1][2].text(-0.49, 0.9,
+  "Vacuum test with $\\gamma={0:.3f}$ in 1D at $t = {1:.2f}$".format(
+    gamma, time), fontsize = 10)
+ax[1][2].text(-0.49, 0.8,
+  "Left:~~ $(P_L, \\rho_L, v_L) = ({0:.3f}, {1:.3f}, {2:.3f})$".format(
+    PL, rhoL, vL), fontsize = 10)
+ax[1][2].text(-0.49, 0.7,
+  "Right: $(P_R, \\rho_R, v_R) = ({0:.3f}, {1:.3f}, {2:.3f})$".format(
+    PR, rhoR, vR), fontsize = 10)
+ax[1][2].plot([-0.49, 0.1], [0.62, 0.62], "k-", lw = 1)
+ax[1][2].text(-0.49, 0.5, "$\\textsc{{Swift}}$ {0}".format(git), fontsize = 10)
+ax[1][2].text(-0.49, 0.4, scheme, fontsize = 10)
+ax[1][2].text(-0.49, 0.3, kernel, fontsize = 10)
+ax[1][2].text(-0.49, 0.2,
+  "${0:.2f}$ neighbours ($\\eta={1:.3f}$)".format(neighbours, eta),
+  fontsize = 10)
+ax[1][2].set_xlim(-0.5, 0.5)
+ax[1][2].set_ylim(0., 1.)
+ax[1][2].set_xticks([])
+ax[1][2].set_yticks([])
+
+pl.tight_layout()
+pl.savefig("Vacuum.png", dpi = 200)

examples/VacuumSpherical_2D/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e glassPlane_128.hdf5 ]
+then
+    echo "Fetching initial glass file for the 2D vacuum expansion example..."
+    ./getGlass.sh
+fi
+if [ ! -e vacuum.hdf5 ]
+then
+    echo "Generating initial conditions for the 2D vacuum expansion example..."
+    python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 4 vacuum.yml 2>&1 | tee output.log
+
+# Get the 1D high resolution reference result if not present.
+if [ ! -e vacuumSpherical2D_exact.txt ]
+then
+    echo "Fetching reference solution for the 2D vacuum expansion test..."
+    ./getReference.sh
+fi
+
+# Plot the result
+python plotSolution.py 1

examples/VacuumSpherical_2D/vacuum.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.05   # The end time of the simulation (in internal units).
+  dt_min:     1e-10  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            vacuum # Common part of the name of output files
+  time_first:          0.       # Time of the first output (in internal units)
+  delta_time:          0.05     # Time difference between consecutive outputs (in internal units)
+
+# 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.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./vacuum.hdf5       # The file to read
+

examples/VacuumSpherical_3D/getGlass.sh

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

examples/VacuumSpherical_3D/getReference.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/vacuumSpherical3D_exact.txt

examples/VacuumSpherical_3D/makeIC.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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 numpy as np
+import h5py
+
+# Generates an overdensity within a vacuum to test the vacuum resolving
+# capabilities of the code
+
+# Parameters
+gamma = 5. / 3.    # Gas adiabatic index
+
+fileName = "vacuum.hdf5" 
+
+#---------------------------------------------------
+glass = h5py.File("glassCube_64.hdf5", "r")
+
+# Read particle positions and h from the glass
+pos = glass["/PartType0/Coordinates"][:,:]
+h = glass["/PartType0/SmoothingLength"][:] * 0.3
+
+# Make 8 copies of the glass to get more particles
+pos *= 0.5
+h *= 0.5
+pos = np.append(pos, pos + np.array([0.5, 0., 0.]), axis = 0)
+pos = np.append(pos, pos + np.array([0., 0.5, 0.]), axis = 0)
+pos = np.append(pos, pos + np.array([0., 0., 0.5]), axis = 0)
+h = np.append(h, h)
+h = np.append(h, h)
+h = np.append(h, h)
+
+radius = np.sqrt((pos[:,0] - 0.5)**2 + (pos[:,1] - 0.5)**2 + \
+                 (pos[:,2] - 0.5)**2)
+index = radius < 0.25
+pos = pos[index]
+h = h[index]
+
+numPart = len(h)
+vol = 4. * np.pi / 3. * 0.25**3
+
+# Generate extra arrays
+v = np.zeros((numPart, 3))
+ids = np.linspace(1, numPart, numPart)
+m = np.zeros(numPart)
+u = np.zeros(numPart)
+
+m[:] = 1. * vol / numPart
+u[:] = 1. / (1. * (gamma - 1.))
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [1., 1., 1.]
+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
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#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/VacuumSpherical_3D/plotSolution.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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 numpy as np
+import matplotlib
+matplotlib.use("Agg")
+import pylab as pl
+import h5py
+import sys
+import scipy.stats as stats
+
+# Parameters
+gamma = 5. / 3. # Polytropic index
+rhoL = 1.       # Initial density in the non vacuum state
+vL = 0.         # Initial velocity in the non vacuum state
+PL = 1.         # Initial pressure in the non vacuum state
+rhoR = 0.       # Initial vacuum density
+vR = 0.         # Initial vacuum velocity
+PR = 0.         # Initial vacuum pressure
+
+# 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
+}
+pl.rcParams.update(params)
+pl.rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+# Read the snapshot index from the command line argument
+snap = int(sys.argv[1])
+
+# Open the file and read the relevant data
+file = h5py.File("vacuum_{0:04d}.hdf5".format(snap), "r")
+coords = file["/PartType0/Coordinates"]
+x = np.sqrt((coords[:,0] - 0.5)**2 + (coords[:,1] - 0.5)**2 + \
+            (coords[:,2] - 0.5)**2)
+rho = file["/PartType0/Density"][:]
+vels = file["/PartType0/Velocities"]
+v = np.sqrt(vels[:,0]**2 + vels[:,1]**2 + vels[:,2]**2)
+u = file["/PartType0/InternalEnergy"][:]
+S = file["/PartType0/Entropy"][:]
+P = file["/PartType0/Pressure"][:]
+time = file["/Header"].attrs["Time"][0]
+
+scheme = file["/HydroScheme"].attrs["Scheme"]
+kernel = file["/HydroScheme"].attrs["Kernel function"]
+neighbours = file["/HydroScheme"].attrs["Kernel target N_ngb"][0]
+eta = file["/HydroScheme"].attrs["Kernel eta"][0]
+git = file["Code"].attrs["Git Revision"]
+
+# Bin the data values
+# We let scipy choose the bins and then reuse them for all other quantities
+rho_bin, x_bin_edge, _ = \
+  stats.binned_statistic(x, rho, statistic = "mean", bins = 50)
+rho2_bin, _, _ = \
+  stats.binned_statistic(x, rho**2, statistic = "mean", bins = x_bin_edge)
+rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
+
+v_bin, _, _ = \
+  stats.binned_statistic(x, v, statistic = "mean", bins = x_bin_edge)
+v2_bin, _, _ = \
+  stats.binned_statistic(x, v**2, statistic = "mean", bins = x_bin_edge)
+v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
+
+P_bin, _, _ = \
+  stats.binned_statistic(x, P, statistic = "mean", bins = x_bin_edge)
+P2_bin, _, _ = \
+  stats.binned_statistic(x, P**2, statistic = "mean", bins = x_bin_edge)
+P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
+
+u_bin, _, _ = \
+  stats.binned_statistic(x, u, statistic = "mean", bins = x_bin_edge)
+u2_bin, _, _ = \
+  stats.binned_statistic(x, u**2, statistic = "mean", bins = x_bin_edge)
+u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
+
+S_bin, _, _ = \
+  stats.binned_statistic(x, S, statistic = "mean", bins = x_bin_edge)
+S2_bin, _, _ = \
+  stats.binned_statistic(x, S**2, statistic = "mean", bins = x_bin_edge)
+S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
+
+x_bin = 0.5 * (x_bin_edge[1:] + x_bin_edge[:-1])
+
+ref = np.loadtxt("vacuumSpherical3D_exact.txt")
+
+# Plot the interesting quantities
+fig, ax = pl.subplots(2, 3)
+
+# Velocity profile
+ax[0][0].plot(x, v, "r.", markersize = 0.2)
+ax[0][0].plot(ref[:,0], ref[:,2], "k--", alpha = 0.8, linewidth = 1.2)
+ax[0][0].errorbar(x_bin, v_bin, yerr = v_sigma_bin, fmt = ".",
+  markersize = 8., color = "b", linewidth = 1.2)
+ax[0][0].set_xlabel("${\\rm{Radius}}~r$", labelpad = 0)
+ax[0][0].set_ylabel("${\\rm{Velocity}}~v_r$", labelpad = 0)
+ax[0][0].set_xlim(0., 0.4)
+ax[0][0].set_ylim(-0.1, 3.2)
+
+# Density profile
+ax[0][1].plot(x, rho, "r.", markersize = 0.2)
+ax[0][1].plot(ref[:,0], ref[:,1], "k--", alpha = 0.8, linewidth = 1.2)
+ax[0][1].errorbar(x_bin, rho_bin, yerr = rho_sigma_bin, fmt = ".",
+  markersize = 8., color = "b", linewidth = 1.2)
+ax[0][1].set_xlabel("${\\rm{Radius}}~r$", labelpad = 0)
+ax[0][1].set_ylabel("${\\rm{Density}}~\\rho$", labelpad = 0)
+ax[0][1].set_xlim(0., 0.4)
+
+# Pressure profile
+ax[0][2].plot(x, P, "r.", markersize = 0.2)
+ax[0][2].plot(ref[:,0], ref[:,3], "k--", alpha = 0.8, linewidth = 1.2)
+ax[0][2].errorbar(x_bin, P_bin, yerr = P_sigma_bin, fmt = ".",
+  markersize = 8., color = "b", linewidth = 1.2)
+ax[0][2].set_xlabel("${\\rm{Radius}}~r$", labelpad = 0)
+ax[0][2].set_ylabel("${\\rm{Pressure}}~P$", labelpad = 0)
+ax[0][2].set_xlim(0., 0.4)
+
+# Internal energy profile
+ax[1][0].plot(x, u, "r.", markersize = 0.2)
+ax[1][0].plot(ref[:,0], ref[:,3] / ref[:,1] / (gamma - 1.), "k--", alpha = 0.8,
+              linewidth = 1.2)
+ax[1][0].errorbar(x_bin, u_bin, yerr = u_sigma_bin, fmt = ".",
+  markersize = 8., color = "b", linewidth = 1.2)
+ax[1][0].set_xlabel("${\\rm{Radius}}~r$", labelpad = 0)
+ax[1][0].set_ylabel("${\\rm{Internal~Energy}}~u$", labelpad = 0)
+ax[1][0].set_xlim(0., 0.4)
+
+# Entropy profile
+ax[1][1].plot(x, S, "r.", markersize = 0.2)
+ax[1][1].plot(ref[:,0], ref[:,3] / ref[:,1]**gamma, "k--", alpha = 0.8,
+              linewidth = 1.2)
+ax[1][1].errorbar(x_bin, S_bin, yerr = S_sigma_bin, fmt = ".",
+  markersize = 8., color = "b", linewidth = 1.2)
+ax[1][1].set_xlabel("${\\rm{Radius}}~r$", labelpad = 0)
+ax[1][1].set_ylabel("${\\rm{Entropy}}~S$", labelpad = 0)
+ax[1][1].set_xlim(0., 0.4)
+ax[1][1].set_ylim(0., 4.)
+
+# Run information
+ax[1][2].set_frame_on(False)
+ax[1][2].text(-0.49, 0.9,
+  "Vacuum test with $\\gamma={0:.3f}$ in 1D at $t = {1:.2f}$".format(
+    gamma, time), fontsize = 10)
+ax[1][2].text(-0.49, 0.8,
+  "Left:~~ $(P_L, \\rho_L, v_L) = ({0:.3f}, {1:.3f}, {2:.3f})$".format(
+    PL, rhoL, vL), fontsize = 10)
+ax[1][2].text(-0.49, 0.7,
+  "Right: $(P_R, \\rho_R, v_R) = ({0:.3f}, {1:.3f}, {2:.3f})$".format(
+    PR, rhoR, vR), fontsize = 10)
+ax[1][2].plot([-0.49, 0.1], [0.62, 0.62], "k-", lw = 1)
+ax[1][2].text(-0.49, 0.5, "$\\textsc{{Swift}}$ {0}".format(git), fontsize = 10)
+ax[1][2].text(-0.49, 0.4, scheme, fontsize = 10)
+ax[1][2].text(-0.49, 0.3, kernel, fontsize = 10)
+ax[1][2].text(-0.49, 0.2,
+  "${0:.2f}$ neighbours ($\\eta={1:.3f}$)".format(neighbours, eta),
+  fontsize = 10)
+ax[1][2].set_xlim(-0.5, 0.5)
+ax[1][2].set_ylim(0., 1.)
+ax[1][2].set_xticks([])
+ax[1][2].set_yticks([])
+
+pl.tight_layout()
+pl.savefig("Vacuum.png", dpi = 200)

examples/VacuumSpherical_3D/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e glassCube_64.hdf5 ]
+then
+    echo "Fetching initial glass file for the 3D vacuum expansion example..."
+    ./getGlass.sh
+fi
+if [ ! -e vacuum.hdf5 ]
+then
+    echo "Generating initial conditions for the 3D vacuum expansion example..."
+    python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 16 vacuum.yml 2>&1 | tee output.log
+
+# Get the reference solution if it is not present.
+if [ ! -e vacuumSpherical3D_exact.txt ]
+then
+    echo "Fetching reference solution for the 3D vacuum expansion test..."
+    ./getReference.sh
+fi
+
+# Plot the result
+python plotSolution.py 1

examples/VacuumSpherical_3D/vacuum.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.05   # The end time of the simulation (in internal units).
+  dt_min:     1e-10  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            vacuum # Common part of the name of output files
+  time_first:          0.       # Time of the first output (in internal units)
+  delta_time:          0.05     # Time difference between consecutive outputs (in internal units)
+
+# 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.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./vacuum.hdf5       # The file to read
+

src/const.h

@@ -75,7 +75,7 @@
 /* Show a warning message if a pathological configuration has been detected. */
 //#define GIZMO_PATHOLOGICAL_WARNING
 /* Crash if a pathological configuration has been detected. */
-#define GIZMO_PATHOLOGICAL_ERROR
+//#define GIZMO_PATHOLOGICAL_ERROR
 /* Maximum allowed gradient matrix condition number. If the condition number of
    the gradient matrix (defined in equation C1 in Hopkins, 2015) is larger than
    this value, we artificially increase the number of neighbours to get a more