diff --git a/examples/EvrardCollapse_3D/plotSolution.py b/examples/EvrardCollapse_3D/plotSolution.py
index 328417c1bb7dae871cf15d4c512dc7a2a8cdf38d..8422b9c45fd573f3d0ae36324d6e39ab23cceb25 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 0000000000000000000000000000000000000000..9b6e884850571f7dfbc3f2018235118447bc00be
--- /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 342b67eae1061a9a94675a59f53ebf9e4ba7717f..1719162dec34626d6f4ecb8267c4d06f85b3db26 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 65a1a109857da6e32663451430b3284aa7207d5a..31717bd806ddd6c98c24dfc1def6f79dddff42ff 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 f4e7e46c5315106ea019e81ade27a5b5a9e5b4f3..57b7f7ddc64bc25df031eb0cba7547f40d46b31a 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 4f6ae6cd78c29c03aad7639da4d5ac9715052802..539bfba799e3231bd26ae2eb39c271baa1fa6a4b 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 0000000000000000000000000000000000000000..ae3c977064f5e7a408aa249c5fd9089b3c52ecb1
--- /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 0000000000000000000000000000000000000000..01080400dba430852f0245d42fb5e30f971721a8
--- /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 0000000000000000000000000000000000000000..498f1b5bc5277188d8ff8d34a5ec24cd314332d4
--- /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 0000000000000000000000000000000000000000..6a65206ae20ccf79392054d047ba6be04f169f3e
--- /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 0000000000000000000000000000000000000000..51d32b4de679877741b7ecd74238fecb785579e7
--- /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 0000000000000000000000000000000000000000..881b155b62c7f1f2af12a1d013ff5c05f1c16a88
--- /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 0000000000000000000000000000000000000000..d5c5f590ac37c9c9431d626a2ea61b0c12c1513c
--- /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 0000000000000000000000000000000000000000..49784c313f0c8758d88970169c30d39c58745ec4
--- /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 0000000000000000000000000000000000000000..d67a30707a904268a09641210a6a3bfcbf305dad
--- /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 0000000000000000000000000000000000000000..c73e48ee2d311692cdf4aa3b0e52f4766b339df8
--- /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 0000000000000000000000000000000000000000..a136929678f745f6a3d0859ba146e1bc1c6c43d0
--- /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 0000000000000000000000000000000000000000..881b155b62c7f1f2af12a1d013ff5c05f1c16a88
--- /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 764251704aedc1258b0eaff2b6220a885e5de2e5..85d64f5f3bf97b3f6ca6df8f230e0e0990687d34 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