From aaae33e774867fd90aa16684860eff2115b3f98d Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@gmail.com>
Date: Fri, 23 Feb 2018 10:32:34 +0000
Subject: [PATCH] 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    |   6 +-
 .../InteractingBlastWaves_1D/getReference.sh  |   2 +
 .../InteractingBlastWaves_1D/plotSolution.py  | 115 ++++++++++-
 examples/InteractingBlastWaves_1D/run.sh      |   9 +-
 examples/SodShockSpherical_2D/plotSolution.py |   2 +-
 examples/SodShockSpherical_3D/plotSolution.py |   2 +-
 examples/VacuumSpherical_2D/getGlass.sh       |   2 +
 examples/VacuumSpherical_2D/getReference.sh   |   2 +
 examples/VacuumSpherical_2D/makeIC.py         | 101 +++++++++
 examples/VacuumSpherical_2D/plotSolution.py   | 192 +++++++++++++++++
 examples/VacuumSpherical_2D/run.sh            |  26 +++
 examples/VacuumSpherical_2D/vacuum.yml        |  34 +++
 examples/VacuumSpherical_3D/getGlass.sh       |   2 +
 examples/VacuumSpherical_3D/getReference.sh   |   2 +
 examples/VacuumSpherical_3D/makeIC.py         | 104 ++++++++++
 examples/VacuumSpherical_3D/plotSolution.py   | 193 ++++++++++++++++++
 examples/VacuumSpherical_3D/run.sh            |  26 +++
 examples/VacuumSpherical_3D/vacuum.yml        |  34 +++
 src/const.h                                   |   2 +-
 19 files changed, 839 insertions(+), 17 deletions(-)
 create mode 100755 examples/InteractingBlastWaves_1D/getReference.sh
 create mode 100755 examples/VacuumSpherical_2D/getGlass.sh
 create mode 100755 examples/VacuumSpherical_2D/getReference.sh
 create mode 100644 examples/VacuumSpherical_2D/makeIC.py
 create mode 100644 examples/VacuumSpherical_2D/plotSolution.py
 create mode 100755 examples/VacuumSpherical_2D/run.sh
 create mode 100644 examples/VacuumSpherical_2D/vacuum.yml
 create mode 100755 examples/VacuumSpherical_3D/getGlass.sh
 create mode 100755 examples/VacuumSpherical_3D/getReference.sh
 create mode 100644 examples/VacuumSpherical_3D/makeIC.py
 create mode 100644 examples/VacuumSpherical_3D/plotSolution.py
 create mode 100755 examples/VacuumSpherical_3D/run.sh
 create mode 100644 examples/VacuumSpherical_3D/vacuum.yml

diff --git a/examples/EvrardCollapse_3D/plotSolution.py b/examples/EvrardCollapse_3D/plotSolution.py
index 328417c1bb..8422b9c45f 100644
--- a/examples/EvrardCollapse_3D/plotSolution.py
+++ b/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)
diff --git a/examples/InteractingBlastWaves_1D/getReference.sh b/examples/InteractingBlastWaves_1D/getReference.sh
new file mode 100755
index 0000000000..9b6e884850
--- /dev/null
+++ b/examples/InteractingBlastWaves_1D/getReference.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/interactingBlastWaves1D_exact.txt
diff --git a/examples/InteractingBlastWaves_1D/plotSolution.py b/examples/InteractingBlastWaves_1D/plotSolution.py
index 342b67eae1..1719162dec 100644
--- a/examples/InteractingBlastWaves_1D/plotSolution.py
+++ b/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)
diff --git a/examples/InteractingBlastWaves_1D/run.sh b/examples/InteractingBlastWaves_1D/run.sh
index 65a1a10985..31717bd806 100755
--- a/examples/InteractingBlastWaves_1D/run.sh
+++ b/examples/InteractingBlastWaves_1D/run.sh
@@ -1,6 +1,6 @@
 #!/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
diff --git a/examples/SodShockSpherical_2D/plotSolution.py b/examples/SodShockSpherical_2D/plotSolution.py
index f4e7e46c53..57b7f7ddc6 100644
--- a/examples/SodShockSpherical_2D/plotSolution.py
+++ b/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)
diff --git a/examples/SodShockSpherical_3D/plotSolution.py b/examples/SodShockSpherical_3D/plotSolution.py
index 4f6ae6cd78..539bfba799 100644
--- a/examples/SodShockSpherical_3D/plotSolution.py
+++ b/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)
diff --git a/examples/VacuumSpherical_2D/getGlass.sh b/examples/VacuumSpherical_2D/getGlass.sh
new file mode 100755
index 0000000000..ae3c977064
--- /dev/null
+++ b/examples/VacuumSpherical_2D/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_128.hdf5
diff --git a/examples/VacuumSpherical_2D/getReference.sh b/examples/VacuumSpherical_2D/getReference.sh
new file mode 100755
index 0000000000..01080400db
--- /dev/null
+++ b/examples/VacuumSpherical_2D/getReference.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/vacuumSpherical2D_exact.txt
diff --git a/examples/VacuumSpherical_2D/makeIC.py b/examples/VacuumSpherical_2D/makeIC.py
new file mode 100644
index 0000000000..498f1b5bc5
--- /dev/null
+++ b/examples/VacuumSpherical_2D/makeIC.py
@@ -0,0 +1,101 @@
+###############################################################################
+# 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()
diff --git a/examples/VacuumSpherical_2D/plotSolution.py b/examples/VacuumSpherical_2D/plotSolution.py
new file mode 100644
index 0000000000..6a65206ae2
--- /dev/null
+++ b/examples/VacuumSpherical_2D/plotSolution.py
@@ -0,0 +1,192 @@
+###############################################################################
+# 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)
diff --git a/examples/VacuumSpherical_2D/run.sh b/examples/VacuumSpherical_2D/run.sh
new file mode 100755
index 0000000000..51d32b4de6
--- /dev/null
+++ b/examples/VacuumSpherical_2D/run.sh
@@ -0,0 +1,26 @@
+#!/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
diff --git a/examples/VacuumSpherical_2D/vacuum.yml b/examples/VacuumSpherical_2D/vacuum.yml
new file mode 100644
index 0000000000..881b155b62
--- /dev/null
+++ b/examples/VacuumSpherical_2D/vacuum.yml
@@ -0,0 +1,34 @@
+# 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
+
diff --git a/examples/VacuumSpherical_3D/getGlass.sh b/examples/VacuumSpherical_3D/getGlass.sh
new file mode 100755
index 0000000000..d5c5f590ac
--- /dev/null
+++ b/examples/VacuumSpherical_3D/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
diff --git a/examples/VacuumSpherical_3D/getReference.sh b/examples/VacuumSpherical_3D/getReference.sh
new file mode 100755
index 0000000000..49784c313f
--- /dev/null
+++ b/examples/VacuumSpherical_3D/getReference.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/vacuumSpherical3D_exact.txt
diff --git a/examples/VacuumSpherical_3D/makeIC.py b/examples/VacuumSpherical_3D/makeIC.py
new file mode 100644
index 0000000000..d67a30707a
--- /dev/null
+++ b/examples/VacuumSpherical_3D/makeIC.py
@@ -0,0 +1,104 @@
+###############################################################################
+# 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()
diff --git a/examples/VacuumSpherical_3D/plotSolution.py b/examples/VacuumSpherical_3D/plotSolution.py
new file mode 100644
index 0000000000..c73e48ee2d
--- /dev/null
+++ b/examples/VacuumSpherical_3D/plotSolution.py
@@ -0,0 +1,193 @@
+###############################################################################
+# 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)
diff --git a/examples/VacuumSpherical_3D/run.sh b/examples/VacuumSpherical_3D/run.sh
new file mode 100755
index 0000000000..a136929678
--- /dev/null
+++ b/examples/VacuumSpherical_3D/run.sh
@@ -0,0 +1,26 @@
+#!/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
diff --git a/examples/VacuumSpherical_3D/vacuum.yml b/examples/VacuumSpherical_3D/vacuum.yml
new file mode 100644
index 0000000000..881b155b62
--- /dev/null
+++ b/examples/VacuumSpherical_3D/vacuum.yml
@@ -0,0 +1,34 @@
+# 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
+
diff --git a/src/const.h b/src/const.h
index 764251704a..85d64f5f3b 100644
--- a/src/const.h
+++ b/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
-- 
GitLab