diff --git a/examples/ComovingSodShock_1D/makeIC.py b/examples/ComovingSodShock_1D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..371acc685cd00bd211e024ecc17c976ea5a08f68
--- /dev/null
+++ b/examples/ComovingSodShock_1D/makeIC.py
@@ -0,0 +1,123 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 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 h5py
+from numpy import *
+
+# Generates a swift IC file for the 1D Sod Shock in a periodic box
+
+unit_l_in_cgs = 3.086e18
+unit_m_in_cgs = 2.94e55
+unit_t_in_cgs = 3.086e18
+
+# Parameters
+gamma = 5./3. # Gas adiabatic index
+numPart_L = 800 # Number of particles in the left state
+x_min = -1.
+x_max = 1.
+rho_L = 1. # Density left state
+rho_R = 0.125 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
+a_beg = 0.001
+fileName = "sodShock.hdf5"
+
+
+#---------------------------------------------------
+
+# Find how many particles we actually have
+boxSize = x_max - x_min
+numPart_R = int(numPart_L * (rho_R / rho_L))
+numPart = numPart_L + numPart_R
+
+# Now get the distances
+delta_L = (boxSize/2) / numPart_L
+delta_R = (boxSize/2) / numPart_R
+offset_L = delta_L / 2
+offset_R = delta_R / 2
+
+# Build the arrays
+coords = zeros((numPart, 3))
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+h = zeros(numPart)
+u = zeros(numPart)
+
+# Set the particles on the left
+for i in range(numPart_L):
+ coords[i,0] = x_min + offset_L + i * delta_L
+ u[i] = P_L / (rho_L * (gamma - 1.))
+ h[i] = 1.2348 * delta_L
+ m[i] = boxSize * rho_L / (2. * numPart_L)
+ v[i,0] = v_L
+
+# Set the particles on the right
+for j in range(numPart_R):
+ i = numPart_L + j
+ coords[i,0] = offset_R + j * delta_R
+ u[i] = P_R / (rho_R * (gamma - 1.))
+ h[i] = 1.2348 * delta_R
+ m[i] = boxSize * rho_R / (2. * numPart_R)
+ v[i,0] = v_R
+
+# Shift particles
+coords[:,0] -= x_min
+
+u /= (a_beg**(3. * (gamma - 1.)))
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+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"] = 1
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = unit_l_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = unit_m_in_cgs
+grp.attrs["Unit time in cgs (U_t)"] = unit_t_in_cgs
+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=coords, 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/ComovingSodShock_1D/plotSolution.py b/examples/ComovingSodShock_1D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..95674c04bfafd0cd549b69814df82f9a4f80a949
--- /dev/null
+++ b/examples/ComovingSodShock_1D/plotSolution.py
@@ -0,0 +1,310 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 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/>.
+ #
+ ##############################################################################
+
+# Computes the analytical solution of the Sod shock and plots the SPH answer
+
+
+# Generates the analytical solution for the Sod shock test case
+# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
+# This follows the solution given in (Toro, 2009)
+
+
+# Parameters
+gas_gamma = 5./3. # Polytropic index
+rho_L = 1. # Density left state
+rho_R = 0.125 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
+
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left' : 0.045,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.05,
+'figure.subplot.top' : 0.99,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+snap = int(sys.argv[1])
+
+
+# Read the simulation data
+sim = h5py.File("sodShock_%04d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+anow = sim["/Header"].attrs["Scale-factor"]
+a_i = sim["/Cosmology"].attrs["a_beg"]
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"]
+time = 2. * (1. / np.sqrt(a_i) - 1. / np.sqrt(anow)) / H_0
+scheme = str(sim["/HydroScheme"].attrs["Scheme"])
+kernel = str(sim["/HydroScheme"].attrs["Kernel function"])
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = str(sim["Code"].attrs["Git Revision"])
+
+x = sim["/PartType0/Coordinates"][:,0]
+v = sim["/PartType0/Velocities"][:,0] * anow
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Density"][:]
+try:
+ alpha = sim["/PartType0/Viscosity"][:]
+ plot_alpha = True
+except:
+ plot_alpha = False
+
+N = 1000 # Number of points
+x_min = -1.
+x_max = 1.
+
+x += x_min
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+c_L = sqrt(gas_gamma * P_L / rho_L) # Speed of the rarefaction wave
+c_R = sqrt(gas_gamma * P_R / rho_R) # Speed of the shock front
+
+# Helpful variable
+Gama = (gas_gamma - 1.) / (gas_gamma + 1.)
+beta = (gas_gamma - 1.) / (2. * gas_gamma)
+
+# Characteristic function and its derivative, following Toro (2009)
+def compute_f(P_3, P, c):
+ u = P_3 / P
+ if u > 1:
+ term1 = gas_gamma*((gas_gamma+1.)*u + gas_gamma-1.)
+ term2 = sqrt(2./term1)
+ fp = (u - 1.)*c*term2
+ dfdp = c*term2/P + (u - 1.)*c/term2*(-1./term1**2)*gas_gamma*(gas_gamma+1.)/P
+ else:
+ fp = (u**beta - 1.)*(2.*c/(gas_gamma-1.))
+ dfdp = 2.*c/(gas_gamma-1.)*beta*u**(beta-1.)/P
+ return (fp, dfdp)
+
+# Solution of the Riemann problem following Toro (2009)
+def RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R):
+ P_new = ((c_L + c_R + (v_L - v_R)*0.5*(gas_gamma-1.))/(c_L / P_L**beta + c_R / P_R**beta))**(1./beta)
+ P_3 = 0.5*(P_R + P_L)
+ f_L = 1.
+ while fabs(P_3 - P_new) > 1e-6:
+ P_3 = P_new
+ (f_L, dfdp_L) = compute_f(P_3, P_L, c_L)
+ (f_R, dfdp_R) = compute_f(P_3, P_R, c_R)
+ f = f_L + f_R + (v_R - v_L)
+ df = dfdp_L + dfdp_R
+ dp = -f/df
+ prnew = P_3 + dp
+ v_3 = v_L - f_L
+ return (P_new, v_3)
+
+
+# Solve Riemann problem for post-shock region
+(P_3, v_3) = RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R)
+
+# Check direction of shocks and wave
+shock_R = (P_3 > P_R)
+shock_L = (P_3 > P_L)
+
+# Velocity of shock front and and rarefaction wave
+if shock_R:
+ v_right = v_R + c_R**2*(P_3/P_R - 1.)/(gas_gamma*(v_3-v_R))
+else:
+ v_right = c_R + 0.5*(gas_gamma+1.)*v_3 - 0.5*(gas_gamma-1.)*v_R
+
+if shock_L:
+ v_left = v_L + c_L**2*(P_3/p_L - 1.)/(gas_gamma*(v_3-v_L))
+else:
+ v_left = c_L - 0.5*(gas_gamma+1.)*v_3 + 0.5*(gas_gamma-1.)*v_L
+
+# Compute position of the transitions
+x_23 = -fabs(v_left) * time
+if shock_L :
+ x_12 = -fabs(v_left) * time
+else:
+ x_12 = -(c_L - v_L) * time
+
+x_34 = v_3 * time
+
+x_45 = fabs(v_right) * time
+if shock_R:
+ x_56 = fabs(v_right) * time
+else:
+ x_56 = (c_R + v_R) * time
+
+
+# Prepare arrays
+delta_x = (x_max - x_min) / N
+x_s = arange(x_min, x_max, delta_x)
+rho_s = zeros(N)
+P_s = zeros(N)
+v_s = zeros(N)
+
+# Compute solution in the different regions
+for i in range(N):
+ if x_s[i] <= x_12:
+ rho_s[i] = rho_L
+ P_s[i] = P_L
+ v_s[i] = v_L
+ if x_s[i] >= x_12 and x_s[i] < x_23:
+ if shock_L:
+ rho_s[i] = rho_L*(Gama + P_3/P_L)/(1. + Gama * P_3/P_L)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ else:
+ rho_s[i] = rho_L*(Gama * (0. - x_s[i])/(c_L * time) + Gama * v_L/c_L + (1.-Gama))**(2./(gas_gamma-1.))
+ P_s[i] = P_L*(rho_s[i] / rho_L)**gas_gamma
+ v_s[i] = (1.-Gama)*(c_L -(0. - x_s[i]) / time) + Gama*v_L
+ if x_s[i] >= x_23 and x_s[i] < x_34:
+ if shock_L:
+ rho_s[i] = rho_L*(Gama + P_3/P_L)/(1+Gama * P_3/p_L)
+ else:
+ rho_s[i] = rho_L*(P_3 / P_L)**(1./gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_34 and x_s[i] < x_45:
+ if shock_R:
+ rho_s[i] = rho_R*(Gama + P_3/P_R)/(1. + Gama * P_3/P_R)
+ else:
+ rho_s[i] = rho_R*(P_3 / P_R)**(1./gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_45 and x_s[i] < x_56:
+ if shock_R:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+ else:
+ rho_s[i] = rho_R*(Gama*(x_s[i])/(c_R*time) - Gama*v_R/c_R + (1.-Gama))**(2./(gas_gamma-1.))
+ P_s[i] = p_R*(rho_s[i]/rho_R)**gas_gamma
+ v_s[i] = (1.-Gama)*(-c_R - (-x_s[i])/time) + Gama*v_R
+ if x_s[i] >= x_56:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+
+
+# Additional arrays
+u_s = P_s / (rho_s * (gas_gamma - 1.)) #internal energy
+s_s = P_s / rho_s**gas_gamma # entropic function
+
+
+# Plot the interesting quantities
+figure()
+
+# Velocity profile --------------------------------
+subplot(231)
+plot(x, v, '.', color='r', ms=4.0)
+plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(-0.1, 0.95)
+
+# Density profile --------------------------------
+subplot(232)
+plot(x, rho, '.', color='r', ms=4.0)
+plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.05, 1.1)
+
+# Pressure profile --------------------------------
+subplot(233)
+plot(x, P, '.', color='r', ms=4.0)
+plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.01, 1.1)
+
+# Internal energy profile -------------------------
+subplot(234)
+plot(x, u, '.', color='r', ms=4.0)
+plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 2.2)
+
+# Entropy/alpha profile ---------------------------------
+subplot(235)
+
+if plot_alpha:
+ plot(x, alpha, '.', color='r', ms=4.0)
+ ylabel(r"${\rm{Viscosity}}~\alpha$", labelpad=0)
+ # Show location of shock
+ plot([x_56, x_56], [-100, 100], color="k", alpha=0.5, ls="dashed", lw=1.2)
+ ylim(0, 1)
+else:
+ plot(x, S, '.', color='r', ms=4.0)
+ plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+ ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+ ylim(0.8, 3.8)
+
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+xlim(-0.5, 0.5)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+z_now = 1. / anow - 1.
+text(-0.49, 0.9, "Sod shock with $\\gamma=%.3f$ in 1D at $z=%.2f$"%(gas_gamma,z_now), fontsize=10)
+text(-0.49, 0.8, "Left:~~ $(P_L, \\rho_L, v_L) = (%.3f, %.3f, %.3f)$"%(P_L, rho_L, v_L), fontsize=10)
+text(-0.49, 0.7, "Right: $(P_R, \\rho_R, v_R) = (%.3f, %.3f, %.3f)$"%(P_R, rho_R, v_R), fontsize=10)
+z_i = 1. / a_i - 1.
+text(-0.49, 0.6, "Initial redshift: $%.2f$"%z_i, fontsize=10)
+plot([-0.49, 0.1], [0.52, 0.52], 'k-', lw=1)
+text(-0.49, 0.4, "$\\textsc{Swift}$ %s"%git, fontsize=10)
+text(-0.49, 0.3, scheme, fontsize=10)
+text(-0.49, 0.2, kernel, fontsize=10)
+text(-0.49, 0.1, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+tight_layout()
+
+savefig("SodShock.png", dpi=200)
diff --git a/examples/ComovingSodShock_1D/run.sh b/examples/ComovingSodShock_1D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..0d1fc4f1be8699929b2cf3b2ea2c8813ebef9f10
--- /dev/null
+++ b/examples/ComovingSodShock_1D/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e sodShock.hdf5 ]
+then
+ echo "Generating initial conditions for the 1D SodShock example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift --cosmology --hydro --threads=1 sodShock.yml 2>&1 | tee output.log
+
+# Plot the result
+python plotSolution.py 1
diff --git a/examples/ComovingSodShock_1D/sodShock.yml b/examples/ComovingSodShock_1D/sodShock.yml
new file mode 100644
index 0000000000000000000000000000000000000000..2d7a5727cbbc2cd417527ce05d7a8ea8ea05dd71
--- /dev/null
+++ b/examples/ComovingSodShock_1D/sodShock.yml
@@ -0,0 +1,43 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 2.94e55 # Grams
+ UnitLength_in_cgs: 3.086e18 # pc
+ UnitVelocity_in_cgs: 1. # km per s
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sodShock # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1.06638 # Time difference between consecutive outputs (in internal units)
+ scale_factor_first: 0.001
+ compression: 1
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1.02 # 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: ./sodShock.hdf5 # The file to read
+ periodic: 1
+
+Cosmology:
+ Omega_m: 1.
+ Omega_lambda: 0.
+ Omega_b: 1.
+ h: 1.
+ a_begin: 0.001
+ a_end: 0.00106638
+
diff --git a/examples/ComovingSodShock_2D/getGlass.sh b/examples/ComovingSodShock_2D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..f4cb4ebcb4b452b2b123462bc97eed532f43ba25
--- /dev/null
+++ b/examples/ComovingSodShock_2D/getGlass.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_128.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_48.hdf5
diff --git a/examples/ComovingSodShock_2D/makeIC.py b/examples/ComovingSodShock_2D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..51a408866047534f86fbded071d604ec294ed0b7
--- /dev/null
+++ b/examples/ComovingSodShock_2D/makeIC.py
@@ -0,0 +1,127 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 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 h5py
+from numpy import *
+
+# Generates a swift IC file for the 2D Sod Shock in a periodic box
+
+unit_l_in_cgs = 3.086e18
+unit_m_in_cgs = 2.94e55
+unit_t_in_cgs = 3.086e18
+
+# Parameters
+gamma = 5./3. # Gas adiabatic index
+x_min = -1.
+x_max = 1.
+rho_L = 1. # Density left state
+rho_R = 0.140625 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
+a_beg = 0.001
+fileName = "sodShock.hdf5"
+
+
+#---------------------------------------------------
+boxSize = (x_max - x_min)
+
+glass_L = h5py.File("glassPlane_128.hdf5", "r")
+glass_R = h5py.File("glassPlane_48.hdf5", "r")
+
+pos_L = glass_L["/PartType0/Coordinates"][:,:] * 0.5
+pos_R = glass_R["/PartType0/Coordinates"][:,:] * 0.5
+h_L = glass_L["/PartType0/SmoothingLength"][:] * 0.5
+h_R = glass_R["/PartType0/SmoothingLength"][:] * 0.5
+
+# Merge things
+aa = pos_L - array([0.5, 0., 0.])
+pos_LL = append(pos_L, pos_L + array([0.5, 0., 0.]), axis=0)
+pos_RR = append(pos_R, pos_R + array([0.5, 0., 0.]), axis=0)
+pos = append(pos_LL - array([1.0, 0., 0.]), pos_RR, axis=0)
+h_LL = append(h_L, h_L)
+h_RR = append(h_R, h_R)
+h = append(h_LL, h_RR)
+
+numPart_L = size(h_LL)
+numPart_R = size(h_RR)
+numPart = size(h)
+
+vol_L = 0.5
+vol_R = 0.5
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+u = zeros(numPart)
+
+for i in range(numPart):
+ x = pos[i,0]
+
+ if x < 0: #left
+ u[i] = P_L / (rho_L * (gamma - 1.))
+ m[i] = rho_L * vol_L / numPart_L
+ v[i,0] = v_L
+ else: #right
+ u[i] = P_R / (rho_R * (gamma - 1.))
+ m[i] = rho_R * vol_R / numPart_R
+ v[i,0] = v_R
+
+# Shift particles
+pos[:,0] -= x_min
+
+u /= (a_beg**(3. * (gamma - 1.)))
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [boxSize, 0.5, 1.0]
+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
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = unit_l_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = unit_m_in_cgs
+grp.attrs["Unit time in cgs (U_t)"] = unit_t_in_cgs
+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/ComovingSodShock_2D/plotSolution.py b/examples/ComovingSodShock_2D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..8adb3cf5c550ab9724f6a8f34c1a1260a25712e1
--- /dev/null
+++ b/examples/ComovingSodShock_2D/plotSolution.py
@@ -0,0 +1,318 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 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/>.
+ #
+ ##############################################################################
+
+# Computes the analytical solution of the Sod shock and plots the SPH answer
+
+
+# Generates the analytical solution for the Sod shock test case
+# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
+# This follows the solution given in (Toro, 2009)
+
+
+# Parameters
+gas_gamma = 5./3. # Polytropic index
+rho_L = 1. # Density left state
+rho_R = 0.140625 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
+
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+from scipy import stats
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left' : 0.045,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.05,
+'figure.subplot.top' : 0.99,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+snap = int(sys.argv[1])
+
+
+# Read the simulation data
+sim = h5py.File("sodShock_%04d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+anow = sim["/Header"].attrs["Scale-factor"]
+a_i = sim["/Cosmology"].attrs["a_beg"]
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"]
+time = 2. * (1. / np.sqrt(a_i) - 1. / np.sqrt(anow)) / H_0
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+
+x = sim["/PartType0/Coordinates"][:,0]
+v = sim["/PartType0/Velocities"][:,0] * anow
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Density"][:]
+
+N = 1000 # Number of points
+x_min = -1.
+x_max = 1.
+x += x_min
+
+
+# Bin te data
+x_bin_edge = np.arange(-0.6, 0.6, 0.02)
+x_bin = 0.5*(x_bin_edge[1:] + x_bin_edge[:-1])
+rho_bin,_,_ = stats.binned_statistic(x, rho, statistic='mean', bins=x_bin_edge)
+v_bin,_,_ = stats.binned_statistic(x, v, statistic='mean', bins=x_bin_edge)
+P_bin,_,_ = stats.binned_statistic(x, P, statistic='mean', bins=x_bin_edge)
+S_bin,_,_ = stats.binned_statistic(x, S, statistic='mean', bins=x_bin_edge)
+u_bin,_,_ = stats.binned_statistic(x, u, statistic='mean', bins=x_bin_edge)
+rho2_bin,_,_ = stats.binned_statistic(x, rho**2, statistic='mean', bins=x_bin_edge)
+v2_bin,_,_ = stats.binned_statistic(x, v**2, statistic='mean', bins=x_bin_edge)
+P2_bin,_,_ = stats.binned_statistic(x, P**2, statistic='mean', bins=x_bin_edge)
+S2_bin,_,_ = stats.binned_statistic(x, S**2, statistic='mean', bins=x_bin_edge)
+u2_bin,_,_ = stats.binned_statistic(x, u**2, statistic='mean', bins=x_bin_edge)
+rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
+v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
+P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
+S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
+u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
+
+
+# Analytic solution
+c_L = sqrt(gas_gamma * P_L / rho_L) # Speed of the rarefaction wave
+c_R = sqrt(gas_gamma * P_R / rho_R) # Speed of the shock front
+
+# Helpful variable
+Gama = (gas_gamma - 1.) / (gas_gamma + 1.)
+beta = (gas_gamma - 1.) / (2. * gas_gamma)
+
+# Characteristic function and its derivative, following Toro (2009)
+def compute_f(P_3, P, c):
+ u = P_3 / P
+ if u > 1:
+ term1 = gas_gamma*((gas_gamma+1.)*u + gas_gamma-1.)
+ term2 = sqrt(2./term1)
+ fp = (u - 1.)*c*term2
+ dfdp = c*term2/P + (u - 1.)*c/term2*(-1./term1**2)*gas_gamma*(gas_gamma+1.)/P
+ else:
+ fp = (u**beta - 1.)*(2.*c/(gas_gamma-1.))
+ dfdp = 2.*c/(gas_gamma-1.)*beta*u**(beta-1.)/P
+ return (fp, dfdp)
+
+# Solution of the Riemann problem following Toro (2009)
+def RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R):
+ P_new = ((c_L + c_R + (v_L - v_R)*0.5*(gas_gamma-1.))/(c_L / P_L**beta + c_R / P_R**beta))**(1./beta)
+ P_3 = 0.5*(P_R + P_L)
+ f_L = 1.
+ while fabs(P_3 - P_new) > 1e-6:
+ P_3 = P_new
+ (f_L, dfdp_L) = compute_f(P_3, P_L, c_L)
+ (f_R, dfdp_R) = compute_f(P_3, P_R, c_R)
+ f = f_L + f_R + (v_R - v_L)
+ df = dfdp_L + dfdp_R
+ dp = -f/df
+ prnew = P_3 + dp
+ v_3 = v_L - f_L
+ return (P_new, v_3)
+
+
+# Solve Riemann problem for post-shock region
+(P_3, v_3) = RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R)
+
+# Check direction of shocks and wave
+shock_R = (P_3 > P_R)
+shock_L = (P_3 > P_L)
+
+# Velocity of shock front and and rarefaction wave
+if shock_R:
+ v_right = v_R + c_R**2*(P_3/P_R - 1.)/(gas_gamma*(v_3-v_R))
+else:
+ v_right = c_R + 0.5*(gas_gamma+1.)*v_3 - 0.5*(gas_gamma-1.)*v_R
+
+if shock_L:
+ v_left = v_L + c_L**2*(P_3/p_L - 1.)/(gas_gamma*(v_3-v_L))
+else:
+ v_left = c_L - 0.5*(gas_gamma+1.)*v_3 + 0.5*(gas_gamma-1.)*v_L
+
+# Compute position of the transitions
+x_23 = -fabs(v_left) * time
+if shock_L :
+ x_12 = -fabs(v_left) * time
+else:
+ x_12 = -(c_L - v_L) * time
+
+x_34 = v_3 * time
+
+x_45 = fabs(v_right) * time
+if shock_R:
+ x_56 = fabs(v_right) * time
+else:
+ x_56 = (c_R + v_R) * time
+
+
+# Prepare arrays
+delta_x = (x_max - x_min) / N
+x_s = arange(x_min, x_max, delta_x)
+rho_s = zeros(N)
+P_s = zeros(N)
+v_s = zeros(N)
+
+# Compute solution in the different regions
+for i in range(N):
+ if x_s[i] <= x_12:
+ rho_s[i] = rho_L
+ P_s[i] = P_L
+ v_s[i] = v_L
+ if x_s[i] >= x_12 and x_s[i] < x_23:
+ if shock_L:
+ rho_s[i] = rho_L*(Gama + P_3/P_L)/(1. + Gama * P_3/P_L)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ else:
+ rho_s[i] = rho_L*(Gama * (0. - x_s[i])/(c_L * time) + Gama * v_L/c_L + (1.-Gama))**(2./(gas_gamma-1.))
+ P_s[i] = P_L*(rho_s[i] / rho_L)**gas_gamma
+ v_s[i] = (1.-Gama)*(c_L -(0. - x_s[i]) / time) + Gama*v_L
+ if x_s[i] >= x_23 and x_s[i] < x_34:
+ if shock_L:
+ rho_s[i] = rho_L*(Gama + P_3/P_L)/(1+Gama * P_3/p_L)
+ else:
+ rho_s[i] = rho_L*(P_3 / P_L)**(1./gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_34 and x_s[i] < x_45:
+ if shock_R:
+ rho_s[i] = rho_R*(Gama + P_3/P_R)/(1. + Gama * P_3/P_R)
+ else:
+ rho_s[i] = rho_R*(P_3 / P_R)**(1./gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_45 and x_s[i] < x_56:
+ if shock_R:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+ else:
+ rho_s[i] = rho_R*(Gama*(x_s[i])/(c_R*time) - Gama*v_R/c_R + (1.-Gama))**(2./(gas_gamma-1.))
+ P_s[i] = p_R*(rho_s[i]/rho_R)**gas_gamma
+ v_s[i] = (1.-Gama)*(-c_R - (-x_s[i])/time) + Gama*v_R
+ if x_s[i] >= x_56:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+
+
+# Additional arrays
+u_s = P_s / (rho_s * (gas_gamma - 1.)) #internal energy
+s_s = P_s / rho_s**gas_gamma # entropic function
+
+
+# Plot the interesting quantities
+figure()
+
+# Velocity profile --------------------------------
+subplot(231)
+plot(x, v, '.', color='r', ms=0.2)
+plot(x_s, v_s, '--', color='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{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(-0.1, 0.95)
+
+# Density profile --------------------------------
+subplot(232)
+plot(x, rho, '.', color='r', ms=0.2)
+plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.05, 1.1)
+
+# Pressure profile --------------------------------
+subplot(233)
+plot(x, P, '.', color='r', ms=0.2)
+plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(x_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.01, 1.1)
+
+# Internal energy profile -------------------------
+subplot(234)
+plot(x, u, '.', color='r', ms=0.2)
+plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(x_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 2.2)
+
+# Entropy profile ---------------------------------
+subplot(235)
+plot(x, S, '.', color='r', ms=0.2)
+plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 3.8)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+z_now = 1. / anow - 1.
+text(-0.49, 0.9, "Sod shock with $\\gamma=%.3f$ in 2D at $z=%.2f$"%(gas_gamma,z_now), fontsize=10)
+text(-0.49, 0.8, "Left:~~ $(P_L, \\rho_L, v_L) = (%.3f, %.3f, %.3f)$"%(P_L, rho_L, v_L), fontsize=10)
+text(-0.49, 0.7, "Right: $(P_R, \\rho_R, v_R) = (%.3f, %.3f, %.3f)$"%(P_R, rho_R, v_R), fontsize=10)
+z_i = 1. / a_i - 1.
+text(-0.49, 0.6, "Initial redshift: $%.2f$"%z_i, fontsize=10)
+plot([-0.49, 0.1], [0.52, 0.52], 'k-', lw=1)
+text(-0.49, 0.4, "$\\textsc{Swift}$ %s"%git, fontsize=10)
+text(-0.49, 0.3, scheme, fontsize=10)
+text(-0.49, 0.2, kernel, fontsize=10)
+text(-0.49, 0.1, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+tight_layout()
+savefig("SodShock.png", dpi=200)
diff --git a/examples/ComovingSodShock_2D/run.sh b/examples/ComovingSodShock_2D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..bd6cc317d75a3bd415f074c9eaf48511ab693598
--- /dev/null
+++ b/examples/ComovingSodShock_2D/run.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e glassPlane_128.hdf5 ]
+then
+ echo "Fetching initial glass file for the Sod shock example..."
+ ./getGlass.sh
+fi
+if [ ! -e sodShock.hdf5 ]
+then
+ echo "Generating initial conditions for the Sod shock example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift --cosmology --hydro --threads=4 sodShock.yml 2>&1 | tee output.log
+
+python plotSolution.py 1
diff --git a/examples/ComovingSodShock_2D/sodShock.yml b/examples/ComovingSodShock_2D/sodShock.yml
new file mode 100644
index 0000000000000000000000000000000000000000..2d7a5727cbbc2cd417527ce05d7a8ea8ea05dd71
--- /dev/null
+++ b/examples/ComovingSodShock_2D/sodShock.yml
@@ -0,0 +1,43 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 2.94e55 # Grams
+ UnitLength_in_cgs: 3.086e18 # pc
+ UnitVelocity_in_cgs: 1. # km per s
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sodShock # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1.06638 # Time difference between consecutive outputs (in internal units)
+ scale_factor_first: 0.001
+ compression: 1
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1.02 # 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: ./sodShock.hdf5 # The file to read
+ periodic: 1
+
+Cosmology:
+ Omega_m: 1.
+ Omega_lambda: 0.
+ Omega_b: 1.
+ h: 1.
+ a_begin: 0.001
+ a_end: 0.00106638
+
diff --git a/examples/ComovingSodShock_3D/README.txt b/examples/ComovingSodShock_3D/README.txt
new file mode 100644
index 0000000000000000000000000000000000000000..2b2f0d16207079f5d29e9ec281fe5255cc5886fb
--- /dev/null
+++ b/examples/ComovingSodShock_3D/README.txt
@@ -0,0 +1,20 @@
+Cosmological version of the standard Sod shock test.
+
+In the co-moving coordinates that SWIFT uses, the Euler equations of
+hydrodynamics have an elegant form with a few additional factors that
+involve the scale factor a. For the specific case of a polytropic index
+gamma = 5/3, all additional factors are in fact the same: 1/a^2. For
+this case, hydrodynamics in the co-moving frame is identical to
+hydrodynamics in a physical non-cosmological frame with a rescaled time
+variable dt'=a^2*dt.
+
+We choose an Einstein-de Sitter cosmology with H(a)=H_0*a^(3/2) and a
+box size of 1 pc and rescale the Sod shock initial conditions so that
+the internal coordinates, density and pressure values are still the same
+as for the original Sod shock in the non-cosmological case. We then
+evolve the initial condition from z=999 to z=960.5, which corresponds to
+a rescaled time interval dt'~0.12. If the co-moving coordinates are
+implemented correctly, the resulting co-moving density, internal
+velocity and co-moving pressure profiles should match those for the
+non-co-moving variables in the ordinary non-cosmological Sod shock at
+t=0.12.
diff --git a/examples/ComovingSodShock_3D/getGlass.sh b/examples/ComovingSodShock_3D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..f61b61d4e6c51b44576fd7cdd6242cb9f0133039
--- /dev/null
+++ b/examples/ComovingSodShock_3D/getGlass.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
diff --git a/examples/ComovingSodShock_3D/makeIC.py b/examples/ComovingSodShock_3D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..b6528bc5ab3670d4423945d194fc537c1bb672a1
--- /dev/null
+++ b/examples/ComovingSodShock_3D/makeIC.py
@@ -0,0 +1,126 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 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 h5py
+from numpy import *
+
+# Generates a swift IC file for the 3D Sod Shock in a periodic box
+
+unit_l_in_cgs = 3.086e18
+unit_m_in_cgs = 2.94e55
+unit_t_in_cgs = 3.086e18
+
+# Parameters
+gamma = 5./3. # Gas adiabatic index
+x_min = -1.
+x_max = 1.
+rho_L = 1. # Density left state
+rho_R = 0.125 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
+a_beg = 0.001
+fileName = "sodShock.hdf5"
+
+
+#---------------------------------------------------
+boxSize = (x_max - x_min)
+
+glass_L = h5py.File("glassCube_64.hdf5", "r")
+glass_R = h5py.File("glassCube_32.hdf5", "r")
+
+pos_L = glass_L["/PartType0/Coordinates"][:,:] * 0.5
+pos_R = glass_R["/PartType0/Coordinates"][:,:] * 0.5
+h_L = glass_L["/PartType0/SmoothingLength"][:] * 0.5
+h_R = glass_R["/PartType0/SmoothingLength"][:] * 0.5
+
+# Merge things
+aa = pos_L - array([0.5, 0., 0.])
+pos_LL = append(pos_L, pos_L + array([0.5, 0., 0.]), axis=0)
+pos_RR = append(pos_R, pos_R + array([0.5, 0., 0.]), axis=0)
+pos = append(pos_LL - array([1.0, 0., 0.]), pos_RR, axis=0)
+h_LL = append(h_L, h_L)
+h_RR = append(h_R, h_R)
+h = append(h_LL, h_RR)
+
+numPart_L = size(h_LL)
+numPart_R = size(h_RR)
+numPart = size(h)
+
+vol_L = (0.25 * boxSize)**2 * (0.5 * boxSize)
+vol_R = (0.25 * boxSize)**2 * (0.5 * boxSize)
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+u = zeros(numPart)
+
+for i in range(numPart):
+ x = pos[i,0]
+
+ if x < 0: #left
+ u[i] = P_L / (rho_L * (gamma - 1.))
+ m[i] = rho_L * vol_L / numPart_L
+ v[i,0] = v_L
+ else: #right
+ u[i] = P_R / (rho_R * (gamma - 1.))
+ m[i] = rho_R * vol_R / numPart_R
+ v[i,0] = v_R
+
+# Shift particles
+pos[:,0] -= x_min
+
+u /= (a_beg**(3. * (gamma - 1.)))
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [boxSize, 0.25 * boxSize, 0.25 * boxSize]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 3
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = unit_l_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = unit_m_in_cgs
+grp.attrs["Unit time in cgs (U_t)"] = unit_t_in_cgs
+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/ComovingSodShock_3D/plotSolution.py b/examples/ComovingSodShock_3D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..f05a385e8620b18189d2e7abca8aebb8ae65060e
--- /dev/null
+++ b/examples/ComovingSodShock_3D/plotSolution.py
@@ -0,0 +1,316 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 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/>.
+ #
+ ##############################################################################
+
+# Computes the analytical solution of the Sod shock and plots the SPH answer
+
+
+# Generates the analytical solution for the Sod shock test case
+# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
+# This follows the solution given in (Toro, 2009)
+
+
+# Parameters
+gas_gamma = 5./3. # Polytropic index
+rho_L = 1. # Density left state
+rho_R = 0.125 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
+
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+from scipy import stats
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left' : 0.045,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.05,
+'figure.subplot.top' : 0.99,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+snap = int(sys.argv[1])
+
+
+# Read the simulation data
+sim = h5py.File("sodShock_%04d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+anow = sim["/Header"].attrs["Scale-factor"]
+a_i = sim["/Cosmology"].attrs["a_beg"]
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"]
+time = 2. * (1. / np.sqrt(a_i) - 1. / np.sqrt(anow)) / H_0
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+
+x = sim["/PartType0/Coordinates"][:,0]
+v = sim["/PartType0/Velocities"][:,0] * anow
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Density"][:]
+
+x_min = -1.
+x_max = 1.
+x += x_min
+N = 1000
+
+# Bin te data
+x_bin_edge = np.arange(-0.6, 0.6, 0.02)
+x_bin = 0.5*(x_bin_edge[1:] + x_bin_edge[:-1])
+rho_bin,_,_ = stats.binned_statistic(x, rho, statistic='mean', bins=x_bin_edge)
+v_bin,_,_ = stats.binned_statistic(x, v, statistic='mean', bins=x_bin_edge)
+P_bin,_,_ = stats.binned_statistic(x, P, statistic='mean', bins=x_bin_edge)
+S_bin,_,_ = stats.binned_statistic(x, S, statistic='mean', bins=x_bin_edge)
+u_bin,_,_ = stats.binned_statistic(x, u, statistic='mean', bins=x_bin_edge)
+rho2_bin,_,_ = stats.binned_statistic(x, rho**2, statistic='mean', bins=x_bin_edge)
+v2_bin,_,_ = stats.binned_statistic(x, v**2, statistic='mean', bins=x_bin_edge)
+P2_bin,_,_ = stats.binned_statistic(x, P**2, statistic='mean', bins=x_bin_edge)
+S2_bin,_,_ = stats.binned_statistic(x, S**2, statistic='mean', bins=x_bin_edge)
+u2_bin,_,_ = stats.binned_statistic(x, u**2, statistic='mean', bins=x_bin_edge)
+rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
+v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
+P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
+S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
+u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
+
+
+# Analytic solution
+c_L = sqrt(gas_gamma * P_L / rho_L) # Speed of the rarefaction wave
+c_R = sqrt(gas_gamma * P_R / rho_R) # Speed of the shock front
+
+# Helpful variable
+Gama = (gas_gamma - 1.) / (gas_gamma + 1.)
+beta = (gas_gamma - 1.) / (2. * gas_gamma)
+
+# Characteristic function and its derivative, following Toro (2009)
+def compute_f(P_3, P, c):
+ u = P_3 / P
+ if u > 1:
+ term1 = gas_gamma*((gas_gamma+1.)*u + gas_gamma-1.)
+ term2 = sqrt(2./term1)
+ fp = (u - 1.)*c*term2
+ dfdp = c*term2/P + (u - 1.)*c/term2*(-1./term1**2)*gas_gamma*(gas_gamma+1.)/P
+ else:
+ fp = (u**beta - 1.)*(2.*c/(gas_gamma-1.))
+ dfdp = 2.*c/(gas_gamma-1.)*beta*u**(beta-1.)/P
+ return (fp, dfdp)
+
+# Solution of the Riemann problem following Toro (2009)
+def RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R):
+ P_new = ((c_L + c_R + (v_L - v_R)*0.5*(gas_gamma-1.))/(c_L / P_L**beta + c_R / P_R**beta))**(1./beta)
+ P_3 = 0.5*(P_R + P_L)
+ f_L = 1.
+ while fabs(P_3 - P_new) > 1e-6:
+ P_3 = P_new
+ (f_L, dfdp_L) = compute_f(P_3, P_L, c_L)
+ (f_R, dfdp_R) = compute_f(P_3, P_R, c_R)
+ f = f_L + f_R + (v_R - v_L)
+ df = dfdp_L + dfdp_R
+ dp = -f/df
+ prnew = P_3 + dp
+ v_3 = v_L - f_L
+ return (P_new, v_3)
+
+
+# Solve Riemann problem for post-shock region
+(P_3, v_3) = RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R)
+
+# Check direction of shocks and wave
+shock_R = (P_3 > P_R)
+shock_L = (P_3 > P_L)
+
+# Velocity of shock front and and rarefaction wave
+if shock_R:
+ v_right = v_R + c_R**2*(P_3/P_R - 1.)/(gas_gamma*(v_3-v_R))
+else:
+ v_right = c_R + 0.5*(gas_gamma+1.)*v_3 - 0.5*(gas_gamma-1.)*v_R
+
+if shock_L:
+ v_left = v_L + c_L**2*(P_3/p_L - 1.)/(gas_gamma*(v_3-v_L))
+else:
+ v_left = c_L - 0.5*(gas_gamma+1.)*v_3 + 0.5*(gas_gamma-1.)*v_L
+
+# Compute position of the transitions
+x_23 = -fabs(v_left) * time
+if shock_L :
+ x_12 = -fabs(v_left) * time
+else:
+ x_12 = -(c_L - v_L) * time
+
+x_34 = v_3 * time
+
+x_45 = fabs(v_right) * time
+if shock_R:
+ x_56 = fabs(v_right) * time
+else:
+ x_56 = (c_R + v_R) * time
+
+
+# Prepare arrays
+delta_x = (x_max - x_min) / N
+x_s = arange(x_min, x_max, delta_x)
+rho_s = zeros(N)
+P_s = zeros(N)
+v_s = zeros(N)
+
+# Compute solution in the different regions
+for i in range(N):
+ if x_s[i] <= x_12:
+ rho_s[i] = rho_L
+ P_s[i] = P_L
+ v_s[i] = v_L
+ if x_s[i] >= x_12 and x_s[i] < x_23:
+ if shock_L:
+ rho_s[i] = rho_L*(Gama + P_3/P_L)/(1. + Gama * P_3/P_L)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ else:
+ rho_s[i] = rho_L*(Gama * (0. - x_s[i])/(c_L * time) + Gama * v_L/c_L + (1.-Gama))**(2./(gas_gamma-1.))
+ P_s[i] = P_L*(rho_s[i] / rho_L)**gas_gamma
+ v_s[i] = (1.-Gama)*(c_L -(0. - x_s[i]) / time) + Gama*v_L
+ if x_s[i] >= x_23 and x_s[i] < x_34:
+ if shock_L:
+ rho_s[i] = rho_L*(Gama + P_3/P_L)/(1+Gama * P_3/p_L)
+ else:
+ rho_s[i] = rho_L*(P_3 / P_L)**(1./gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_34 and x_s[i] < x_45:
+ if shock_R:
+ rho_s[i] = rho_R*(Gama + P_3/P_R)/(1. + Gama * P_3/P_R)
+ else:
+ rho_s[i] = rho_R*(P_3 / P_R)**(1./gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_45 and x_s[i] < x_56:
+ if shock_R:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+ else:
+ rho_s[i] = rho_R*(Gama*(x_s[i])/(c_R*time) - Gama*v_R/c_R + (1.-Gama))**(2./(gas_gamma-1.))
+ P_s[i] = p_R*(rho_s[i]/rho_R)**gas_gamma
+ v_s[i] = (1.-Gama)*(-c_R - (-x_s[i])/time) + Gama*v_R
+ if x_s[i] >= x_56:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+
+
+# Additional arrays
+u_s = P_s / (rho_s * (gas_gamma - 1.)) #internal energy
+s_s = P_s / rho_s**gas_gamma # entropic function
+
+# Plot the interesting quantities
+figure()
+
+# Velocity profile --------------------------------
+subplot(231)
+plot(x, v, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, v_s, '--', color='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{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(-0.1, 0.95)
+
+# Density profile --------------------------------
+subplot(232)
+plot(x, rho, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.05, 1.1)
+
+# Pressure profile --------------------------------
+subplot(233)
+plot(x, P, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(x_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.01, 1.1)
+
+# Internal energy profile -------------------------
+subplot(234)
+plot(x, u, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(x_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 2.2)
+
+# Entropy profile ---------------------------------
+subplot(235)
+plot(x, S, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 3.8)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+znow = 1. / anow - 1.
+text(-0.49, 0.9, "Sod shock with $\\gamma=%.3f$ in 3D at $z=%.2f$"%(gas_gamma,znow), fontsize=10)
+text(-0.49, 0.8, "Left:~~ $(P_L, \\rho_L, v_L) = (%.3f, %.3f, %.3f)$"%(P_L, rho_L, v_L), fontsize=10)
+text(-0.49, 0.7, "Right: $(P_R, \\rho_R, v_R) = (%.3f, %.3f, %.3f)$"%(P_R, rho_R, v_R), fontsize=10)
+z_i = 1. / a_i - 1.
+text(-0.49, 0.6, "Initial redshift: $%.2f$"%z_i, fontsize = 10)
+plot([-0.49, 0.1], [0.52, 0.52], 'k-', lw=1)
+text(-0.49, 0.4, "$\\textsc{Swift}$ %s"%git, fontsize=10)
+text(-0.49, 0.3, scheme, fontsize=10)
+text(-0.49, 0.2, kernel, fontsize=10)
+text(-0.49, 0.1, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+tight_layout()
+savefig("SodShock.png", dpi=200)
diff --git a/examples/ComovingSodShock_3D/run.sh b/examples/ComovingSodShock_3D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..00e4f669fdb1347ab1b34cfa11821ca011b73120
--- /dev/null
+++ b/examples/ComovingSodShock_3D/run.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e glassCube_64.hdf5 ]
+then
+ echo "Fetching initial glass file for the Sod shock example..."
+ ./getGlass.sh
+fi
+if [ ! -e sodShock.hdf5 ]
+then
+ echo "Generating initial conditions for the Sod shock example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift --cosmology --hydro --threads=4 sodShock.yml 2>&1 | tee output.log
+
+python plotSolution.py 1
diff --git a/examples/ComovingSodShock_3D/sodShock.yml b/examples/ComovingSodShock_3D/sodShock.yml
new file mode 100644
index 0000000000000000000000000000000000000000..2d7a5727cbbc2cd417527ce05d7a8ea8ea05dd71
--- /dev/null
+++ b/examples/ComovingSodShock_3D/sodShock.yml
@@ -0,0 +1,43 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 2.94e55 # Grams
+ UnitLength_in_cgs: 3.086e18 # pc
+ UnitVelocity_in_cgs: 1. # km per s
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sodShock # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1.06638 # Time difference between consecutive outputs (in internal units)
+ scale_factor_first: 0.001
+ compression: 1
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1.02 # 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: ./sodShock.hdf5 # The file to read
+ periodic: 1
+
+Cosmology:
+ Omega_m: 1.
+ Omega_lambda: 0.
+ Omega_b: 1.
+ h: 1.
+ a_begin: 0.001
+ a_end: 0.00106638
+
diff --git a/src/engine.c b/src/engine.c
index 73b7d9b64e744e08e19145a574a3a2b896ab3995..962ce905e1334537448697c1ccbc4191084c92d7 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -4780,6 +4780,7 @@ void engine_print_policy(struct engine *e) {
* @param e The #engine.
*/
void engine_compute_next_snapshot_time(struct engine *e) {
+
/* Do outputlist file case */
if (e->output_list_snapshots) {
output_list_read_next_time(e->output_list_snapshots, e, "snapshots",
@@ -4800,6 +4801,8 @@ void engine_compute_next_snapshot_time(struct engine *e) {
time = e->a_first_snapshot;
else
time = e->time_first_snapshot;
+
+ int found_snapshot_time = 0;
while (time < time_end) {
/* Output time on the integer timeline */
@@ -4809,7 +4812,10 @@ void engine_compute_next_snapshot_time(struct engine *e) {
e->ti_next_snapshot = (time - e->time_begin) / e->time_base;
/* Found it? */
- if (e->ti_next_snapshot > e->ti_current) break;
+ if (e->ti_next_snapshot > e->ti_current) {
+ found_snapshot_time = 1;
+ break;
+ }
if (e->policy & engine_policy_cosmology)
time *= e->delta_time_snapshot;
@@ -4818,7 +4824,7 @@ void engine_compute_next_snapshot_time(struct engine *e) {
}
/* Deal with last snapshot */
- if (e->ti_next_snapshot >= max_nr_timesteps) {
+ if (!found_snapshot_time) {
e->ti_next_snapshot = -1;
if (e->verbose) message("No further output time.");
} else {
@@ -4864,6 +4870,8 @@ void engine_compute_next_statistics_time(struct engine *e) {
time = e->a_first_statistics;
else
time = e->time_first_statistics;
+
+ int found_stats_time = 0;
while (time < time_end) {
/* Output time on the integer timeline */
@@ -4873,7 +4881,10 @@ void engine_compute_next_statistics_time(struct engine *e) {
e->ti_next_stats = (time - e->time_begin) / e->time_base;
/* Found it? */
- if (e->ti_next_stats > e->ti_current) break;
+ if (e->ti_next_stats > e->ti_current) {
+ found_stats_time = 1;
+ break;
+ }
if (e->policy & engine_policy_cosmology)
time *= e->delta_time_statistics;
@@ -4882,7 +4893,7 @@ void engine_compute_next_statistics_time(struct engine *e) {
}
/* Deal with last statistics */
- if (e->ti_next_stats >= max_nr_timesteps) {
+ if (!found_stats_time) {
e->ti_next_stats = -1;
if (e->verbose) message("No further output time.");
} else {
@@ -4929,6 +4940,8 @@ void engine_compute_next_stf_time(struct engine *e) {
time = e->a_first_stf_output;
else
time = e->time_first_stf_output;
+
+ int found_stf_time = 0;
while (time < time_end) {
/* Output time on the integer timeline */
@@ -4938,7 +4951,10 @@ void engine_compute_next_stf_time(struct engine *e) {
e->ti_next_stf = (time - e->time_begin) / e->time_base;
/* Found it? */
- if (e->ti_next_stf > e->ti_current) break;
+ if (e->ti_next_stf > e->ti_current) {
+ found_stf_time = 1;
+ break;
+ }
if (e->policy & engine_policy_cosmology)
time *= e->delta_time_stf;
@@ -4947,7 +4963,7 @@ void engine_compute_next_stf_time(struct engine *e) {
}
/* Deal with last snapshot */
- if (e->ti_next_stf >= max_nr_timesteps) {
+ if (!found_stf_time) {
e->ti_next_stf = -1;
if (e->verbose) message("No further output time.");
} else {