diff --git a/.gitignore b/.gitignore
index ba0028d163dae64143b0c43b465152e357d042b0..b44b5189190bd6593368c9ce8fbbd19d1d34d6ce 100644
--- a/.gitignore
+++ b/.gitignore
@@ -22,9 +22,6 @@ doc/Doxyfile
examples/swift
examples/swift_mpi
-examples/*.xmf
-examples/used_parameters.yml
-examples/energy.txt
examples/*/*.xmf
examples/*/*.hdf5
examples/*/*.png
diff --git a/examples/GreshoVortex_2D/plotSolution.py b/examples/GreshoVortex_2D/plotSolution.py
index 050bca39a6b7d6f985d630e057a475945471086a..7a86daa6a4e5e1dd80888ceac9a6eb6b08dff443 100644
--- a/examples/GreshoVortex_2D/plotSolution.py
+++ b/examples/GreshoVortex_2D/plotSolution.py
@@ -31,6 +31,7 @@ P0 = 0. # Constant additional pressure (should have no impact on the d
import matplotlib
matplotlib.use("Agg")
from pylab import *
+from scipy import stats
import h5py
# Plot parameters
@@ -104,6 +105,26 @@ u = sim["/PartType0/InternalEnergy"][:]
S = sim["/PartType0/Entropy"][:]
P = sim["/PartType0/Pressure"][:]
+# Bin te data
+r_bin_edge = np.arange(0., 1., 0.02)
+r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
+rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
+v_bin,_,_ = stats.binned_statistic(r, v_phi, statistic='mean', bins=r_bin_edge)
+P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
+S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
+u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
+rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
+v2_bin,_,_ = stats.binned_statistic(r, v_phi**2, statistic='mean', bins=r_bin_edge)
+P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
+S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
+u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_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)
+
+
# Plot the interesting quantities
figure()
@@ -113,6 +134,7 @@ subplot(231)
plot(r, v_phi, '.', color='r', ms=0.5)
plot(solution_r, solution_v_phi, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
@@ -125,6 +147,7 @@ subplot(232)
plot(r, rho, '.', color='r', ms=0.5)
plot(solution_r, solution_rho, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
@@ -138,6 +161,7 @@ subplot(233)
plot(r, P, '.', color='r', ms=0.5)
plot(solution_r, solution_P, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
@@ -150,6 +174,7 @@ subplot(234)
plot(r, u, '.', color='r', ms=0.5)
plot(solution_r, solution_u, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
@@ -163,6 +188,7 @@ subplot(235)
plot(r, S, '.', color='r', ms=0.5)
plot(solution_r, solution_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
diff --git a/examples/Noh_1D/makeIC.py b/examples/Noh_1D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..176f3517455db7a8b0994ac7d1e65fb9cb7419d4
--- /dev/null
+++ b/examples/Noh_1D/makeIC.py
@@ -0,0 +1,90 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+import h5py
+from numpy import *
+
+# Generates a swift IC file for the 1D Noh problem test in a periodic box
+
+# Parameters
+numPart = 1001
+gamma = 5./3. # Gas adiabatic index
+rho0 = 1. # Background density
+P0 = 1.e-6 # Background pressure
+fileName = "noh.hdf5"
+
+#---------------------------------------------------
+coords = zeros((numPart, 3))
+h = zeros(numPart)
+vol = 2.
+
+for i in range(numPart):
+ coords[i,0] = i * vol/numPart + vol/(2.*numPart)
+ h[i] = 1.2348 * vol / numPart
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+u = zeros(numPart)
+
+m[:] = rho0 * vol / numPart
+u[:] = P0 / (rho0 * (gamma - 1))
+v[coords[:,0]<vol/2. ,0] = 1
+v[coords[:,0]>vol/2. ,0] = -1
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [vol, vol, vol]
+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
+
+#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=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/Noh_1D/noh.yml b/examples/Noh_1D/noh.yml
new file mode 100644
index 0000000000000000000000000000000000000000..911ebbdc1ca82f9b5cf2f70841d848008cbc5f6c
--- /dev/null
+++ b/examples/Noh_1D/noh.yml
@@ -0,0 +1,35 @@
+# 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.6 # The end time of the simulation (in internal units).
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: noh # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 5e-2 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-5 # 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).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./noh.hdf5 # The file to read
+
diff --git a/examples/Noh_1D/plotSolution.py b/examples/Noh_1D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..f4916af6e6066d21f76c28b5acef41e1907a83fd
--- /dev/null
+++ b/examples/Noh_1D/plotSolution.py
@@ -0,0 +1,175 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+# Computes the analytical solution of the Noh problem and plots the SPH answer
+
+
+# Parameters
+gas_gamma = 5./3. # Polytropic index
+rho0 = 1. # Background density
+P0 = 1.e-6 # Background pressure
+v0 = 1
+
+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("noh_%03d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+
+x = sim["/PartType0/Coordinates"][:,0]
+v = sim["/PartType0/Velocities"][:,0]
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Density"][:]
+
+N = 1001 # Number of points
+x_min = -1
+x_max = 1
+
+x += x_min
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+x_s = np.arange(0, 2., 2./N) - 1.
+rho_s = np.ones(N) * rho0
+P_s = np.ones(N) * rho0
+v_s = np.ones(N) * v0
+
+# Shock position
+u0 = rho0 * P0 * (gas_gamma-1)
+us = 0.5 * (gas_gamma - 1) * v0
+rs = us * time
+
+# Post-shock values
+rho_s[np.abs(x_s) < rs] = rho0 * (gas_gamma + 1) / (gas_gamma - 1)
+P_s[np.abs(x_s) < rs] = 0.5 * rho0 * v0**2 * (gas_gamma + 1)
+v_s[np.abs(x_s) < rs] = 0.
+
+# Pre-shock values
+rho_s[np.abs(x_s) >= rs] = rho0
+P_s[np.abs(x_s) >= rs] = 0
+v_s[x_s >= rs] = -v0
+v_s[x_s <= -rs] = v0
+
+# 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=-4)
+xlim(-0.5, 0.5)
+ylim(-1.2, 1.2)
+
+# 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.95, 4.4)
+
+# 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.05, 1.8)
+
+# 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.05, 0.8)
+
+# Entropy profile ---------------------------------
+subplot(235)
+plot(x, S, '.', color='r', ms=4.0)
+plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
+xlim(-0.5, 0.5)
+ylim(-0.05, 0.2)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Noh problem with $\\gamma=%.3f$ in 1D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "ICs:~~ $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$"%(1e-6, 1., -1.), fontsize=10)
+plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
+text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
+text(-0.49, 0.4, scheme, fontsize=10)
+text(-0.49, 0.3, kernel, fontsize=10)
+text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+
+savefig("Noh.png", dpi=200)
diff --git a/examples/Noh_1D/run.sh b/examples/Noh_1D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..77788bfa8429e2fbf0502068baa70598acaaa791
--- /dev/null
+++ b/examples/Noh_1D/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e noh.hdf5 ]
+then
+ echo "Generating initial conditions for the Noh problem..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 1 noh.yml 2>&1 | tee output.log
+
+# Plot the solution
+python plotSolution.py 12
diff --git a/examples/Noh_2D/getGlass.sh b/examples/Noh_2D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ae3c977064f5e7a408aa249c5fd9089b3c52ecb1
--- /dev/null
+++ b/examples/Noh_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/Noh_2D/makeIC.py b/examples/Noh_2D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..f7239fa3cd188637df929f86451d20a9978bd1f5
--- /dev/null
+++ b/examples/Noh_2D/makeIC.py
@@ -0,0 +1,98 @@
+
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+import h5py
+from numpy import *
+
+# Generates a swift IC file for the 2D Noh problem in a periodic box
+
+# Parameters
+gamma = 5./3. # Gas adiabatic index
+gamma = 5./3. # Gas adiabatic index
+rho0 = 1. # Background density
+P0 = 1.e-6 # Background pressure
+fileName = "noh.hdf5"
+
+#---------------------------------------------------
+glass = h5py.File("glassPlane_128.hdf5", "r")
+
+vol = 4.
+
+pos = glass["/PartType0/Coordinates"][:,:] * sqrt(vol)
+h = glass["/PartType0/SmoothingLength"][:] * sqrt(vol)
+numPart = size(h)
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+
+m = zeros(numPart)
+u = zeros(numPart)
+m[:] = rho0 * vol / numPart
+u[:] = P0 / (rho0 * (gamma - 1))
+
+# Make radial velocities
+#r = sqrt((pos[:,0]-1.)**2 + (pos[:,1]-1.)**2)
+#theta = arctan2((pos[:,1]-1.), (pos[:,0]-1.))
+v[:,0] = -(pos[:,0] - 1.)
+v[:,1] = -(pos[:,1] - 1.)
+
+norm_v = sqrt(v[:,0]**2 + v[:,1]**2)
+v[:,0] /= norm_v
+v[:,1] /= norm_v
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [sqrt(vol), sqrt(vol), 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/Noh_2D/noh.yml b/examples/Noh_2D/noh.yml
new file mode 100644
index 0000000000000000000000000000000000000000..911ebbdc1ca82f9b5cf2f70841d848008cbc5f6c
--- /dev/null
+++ b/examples/Noh_2D/noh.yml
@@ -0,0 +1,35 @@
+# 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.6 # The end time of the simulation (in internal units).
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: noh # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 5e-2 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-5 # 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).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./noh.hdf5 # The file to read
+
diff --git a/examples/Noh_2D/plotSolution.py b/examples/Noh_2D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..a01a712efd412488aea09c3f3c4e8d68323fc916
--- /dev/null
+++ b/examples/Noh_2D/plotSolution.py
@@ -0,0 +1,198 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+# Computes the analytical solution of the Noh problem and plots the SPH answer
+
+
+# Parameters
+gas_gamma = 5./3. # Polytropic index
+rho0 = 1. # Background density
+P0 = 1.e-6 # Background pressure
+v0 = 1
+
+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("noh_%03d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+
+x = sim["/PartType0/Coordinates"][:,0]
+y = sim["/PartType0/Coordinates"][:,1]
+vx = sim["/PartType0/Velocities"][:,0]
+vy = sim["/PartType0/Velocities"][:,1]
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Density"][:]
+
+r = np.sqrt((x-1)**2 + (y-1)**2)
+v = -np.sqrt(vx**2 + vy**2)
+
+# Bin te data
+r_bin_edge = np.arange(0., 1., 0.02)
+r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
+rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
+v_bin,_,_ = stats.binned_statistic(r, v, statistic='mean', bins=r_bin_edge)
+P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
+S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
+u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
+rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
+v2_bin,_,_ = stats.binned_statistic(r, v**2, statistic='mean', bins=r_bin_edge)
+P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
+S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
+u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_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
+N = 1000 # Number of points
+
+x_s = np.arange(0, 2., 2./N) - 1.
+rho_s = np.ones(N) * rho0
+P_s = np.ones(N) * rho0
+v_s = np.ones(N) * v0
+
+# Shock position
+u0 = rho0 * P0 * (gas_gamma-1)
+us = 0.5 * (gas_gamma - 1) * v0
+rs = us * time
+
+# Post-shock values
+rho_s[np.abs(x_s) < rs] = rho0 * ((gas_gamma + 1) / (gas_gamma - 1))**2
+P_s[np.abs(x_s) < rs] = 0.5 * rho0 * v0**2 * (gas_gamma + 1)**2 / (gas_gamma-1)
+v_s[np.abs(x_s) < rs] = 0.
+
+# Pre-shock values
+rho_s[np.abs(x_s) >= rs] = rho0 * (1 + v0 * time/np.abs(x_s[np.abs(x_s) >=rs]))
+P_s[np.abs(x_s) >= rs] = 0
+v_s[x_s >= rs] = -v0
+v_s[x_s <= -rs] = v0
+
+# 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(r, v, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_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_r$", labelpad=-4)
+xlim(0, 0.5)
+ylim(-1.2, 0.4)
+
+# Density profile --------------------------------
+subplot(232)
+plot(r, rho, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+xlim(0, 0.5)
+ylim(0.95, 19)
+
+# Pressure profile --------------------------------
+subplot(233)
+plot(r, P, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(0, 0.5)
+ylim(-0.5, 11)
+
+# Internal energy profile -------------------------
+subplot(234)
+plot(r, u, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(0, 0.5)
+ylim(-0.05, 0.8)
+
+# Entropy profile ---------------------------------
+subplot(235)
+plot(r, S, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
+xlim(0, 0.5)
+ylim(-0.05, 0.2)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Noh problem with $\\gamma=%.3f$ in 2D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "ICs:~~ $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$"%(1e-6, 1., -1.), fontsize=10)
+plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
+text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
+text(-0.49, 0.4, scheme, fontsize=10)
+text(-0.49, 0.3, kernel, fontsize=10)
+text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+
+savefig("Noh.png", dpi=200)
diff --git a/examples/Noh_2D/run.sh b/examples/Noh_2D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..cff200801018e04ea560bd2c3fbd84057aec2d7c
--- /dev/null
+++ b/examples/Noh_2D/run.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e glassPlane_128.hdf5 ]
+then
+ echo "Fetching initial glass file for the Noh problem..."
+ ./getGlass.sh
+fi
+if [ ! -e noh.hdf5 ]
+then
+ echo "Generating initial conditions for the Noh problem..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 2 noh.yml 2>&1 | tee output.log
+
+# Plot the solution
+python plotSolution.py 12
diff --git a/examples/Noh_3D/getGlass.sh b/examples/Noh_3D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..d5c5f590ac37c9c9431d626a2ea61b0c12c1513c
--- /dev/null
+++ b/examples/Noh_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/Noh_3D/makeIC.py b/examples/Noh_3D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..ec8d46639eecdefd55b917a955f13efc8ffe4126
--- /dev/null
+++ b/examples/Noh_3D/makeIC.py
@@ -0,0 +1,100 @@
+
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+import h5py
+from numpy import *
+
+# Generates a swift IC file for the 3D Noh problem in a periodic box
+
+# Parameters
+gamma = 5./3. # Gas adiabatic index
+gamma = 5./3. # Gas adiabatic index
+rho0 = 1. # Background density
+P0 = 1.e-6 # Background pressure
+fileName = "noh.hdf5"
+
+#---------------------------------------------------
+glass = h5py.File("glassCube_64.hdf5", "r")
+
+vol = 8.
+
+pos = glass["/PartType0/Coordinates"][:,:] * cbrt(vol)
+h = glass["/PartType0/SmoothingLength"][:] * cbrt(vol)
+numPart = size(h)
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+
+m = zeros(numPart)
+u = zeros(numPart)
+m[:] = rho0 * vol / numPart
+u[:] = P0 / (rho0 * (gamma - 1))
+
+# Make radial velocities
+#r = sqrt((pos[:,0]-1.)**2 + (pos[:,1]-1.)**2)
+#theta = arctan2((pos[:,1]-1.), (pos[:,0]-1.))
+v[:,0] = -(pos[:,0] - 1.)
+v[:,1] = -(pos[:,1] - 1.)
+v[:,2] = -(pos[:,2] - 1.)
+
+norm_v = sqrt(v[:,0]**2 + v[:,1]**2 + v[:,2]**2)
+v[:,0] /= norm_v
+v[:,1] /= norm_v
+v[:,2] /= norm_v
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [cbrt(vol), cbrt(vol), cbrt(vol)]
+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/Noh_3D/noh.yml b/examples/Noh_3D/noh.yml
new file mode 100644
index 0000000000000000000000000000000000000000..911ebbdc1ca82f9b5cf2f70841d848008cbc5f6c
--- /dev/null
+++ b/examples/Noh_3D/noh.yml
@@ -0,0 +1,35 @@
+# 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.6 # The end time of the simulation (in internal units).
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: noh # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 5e-2 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-5 # 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).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./noh.hdf5 # The file to read
+
diff --git a/examples/Noh_3D/plotSolution.py b/examples/Noh_3D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..1742e13a5daeff392690a9804fb2831ef4304963
--- /dev/null
+++ b/examples/Noh_3D/plotSolution.py
@@ -0,0 +1,202 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+# Computes the analytical solution of the Noh problem and plots the SPH answer
+
+
+# Parameters
+gas_gamma = 5./3. # Polytropic index
+rho0 = 1. # Background density
+P0 = 1.e-6 # Background pressure
+v0 = 1
+
+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("noh_%03d.hdf5"%snap, "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"]
+kernel = sim["/HydroScheme"].attrs["Kernel function"]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
+eta = sim["/HydroScheme"].attrs["Kernel eta"]
+git = sim["Code"].attrs["Git Revision"]
+
+x = sim["/PartType0/Coordinates"][:,0]
+y = sim["/PartType0/Coordinates"][:,1]
+z = sim["/PartType0/Coordinates"][:,2]
+vx = sim["/PartType0/Velocities"][:,0]
+vy = sim["/PartType0/Velocities"][:,1]
+vz = sim["/PartType0/Velocities"][:,2]
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Density"][:]
+
+r = np.sqrt((x-1)**2 + (y-1)**2 + (z-1)**2)
+v = -np.sqrt(vx**2 + vy**2 + vz**2)
+
+# Bin te data
+r_bin_edge = np.arange(0., 1., 0.02)
+r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
+rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
+v_bin,_,_ = stats.binned_statistic(r, v, statistic='mean', bins=r_bin_edge)
+P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
+S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
+u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
+rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
+v2_bin,_,_ = stats.binned_statistic(r, v**2, statistic='mean', bins=r_bin_edge)
+P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
+S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
+u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_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
+N = 1000 # Number of points
+
+x_s = np.arange(0, 2., 2./N) - 1.
+rho_s = np.ones(N) * rho0
+P_s = np.ones(N) * rho0
+v_s = np.ones(N) * v0
+
+# Shock position
+u0 = rho0 * P0 * (gas_gamma-1)
+us = 0.5 * (gas_gamma - 1) * v0
+rs = us * time
+
+# Post-shock values
+rho_s[np.abs(x_s) < rs] = rho0 * ((gas_gamma + 1) / (gas_gamma - 1))**3
+P_s[np.abs(x_s) < rs] = 0.5 * rho0 * v0**2 * (gas_gamma + 1)**3 / (gas_gamma-1)**2
+v_s[np.abs(x_s) < rs] = 0.
+
+# Pre-shock values
+rho_s[np.abs(x_s) >= rs] = rho0 * (1 + v0 * time/np.abs(x_s[np.abs(x_s) >=rs]))**2
+P_s[np.abs(x_s) >= rs] = 0
+v_s[x_s >= rs] = -v0
+v_s[x_s <= -rs] = v0
+
+# 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(r, v, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_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_r$", labelpad=-4)
+xlim(0, 0.5)
+ylim(-1.2, 0.4)
+
+# Density profile --------------------------------
+subplot(232)
+plot(r, rho, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+xlim(0, 0.5)
+ylim(0.95, 71)
+
+# Pressure profile --------------------------------
+subplot(233)
+plot(r, P, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(0, 0.5)
+ylim(-0.5, 25)
+
+# Internal energy profile -------------------------
+subplot(234)
+plot(r, u, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(0, 0.5)
+ylim(-0.05, 0.8)
+
+# Entropy profile ---------------------------------
+subplot(235)
+plot(r, S, '.', color='r', ms=0.5, alpha=0.2)
+plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Entropy}}~S$", labelpad=-9)
+xlim(0, 0.5)
+ylim(-0.05, 0.2)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Noh problem with $\\gamma=%.3f$ in 3D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "ICs:~~ $(P_0, \\rho_0, v_0) = (%1.2e, %.3f, %.3f)$"%(1e-6, 1., -1.), fontsize=10)
+plot([-0.49, 0.1], [0.62, 0.62], 'k-', lw=1)
+text(-0.49, 0.5, "$\\textsc{Swift}$ %s"%git, fontsize=10)
+text(-0.49, 0.4, scheme, fontsize=10)
+text(-0.49, 0.3, kernel, fontsize=10)
+text(-0.49, 0.2, "$%.2f$ neighbours ($\\eta=%.3f$)"%(neighbours, eta), fontsize=10)
+xlim(-0.5, 0.5)
+ylim(0, 1)
+xticks([])
+yticks([])
+
+
+savefig("Noh.png", dpi=200)
diff --git a/examples/Noh_3D/run.sh b/examples/Noh_3D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..b9e4fb145b2465433aa2bc0362aba19cc1267461
--- /dev/null
+++ b/examples/Noh_3D/run.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e glassCube_64.hdf5 ]
+then
+ echo "Fetching initial glass file for the Noh problem..."
+ ./getGlass.sh
+fi
+if [ ! -e noh.hdf5 ]
+then
+ echo "Generating initial conditions for the Noh problem..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 2 noh.yml 2>&1 | tee output.log
+
+# Plot the solution
+python plotSolution.py 12
diff --git a/examples/SedovBlast_2D/plotSolution.py b/examples/SedovBlast_2D/plotSolution.py
index 69e4e1232dd5c14b06e8a705f4add391f1f597f0..d8c0c9791d1834cc2a5cf0103b46a49e20d2e8a3 100644
--- a/examples/SedovBlast_2D/plotSolution.py
+++ b/examples/SedovBlast_2D/plotSolution.py
@@ -35,6 +35,7 @@ gas_gamma = 5./3. # Gas polytropic index
import matplotlib
matplotlib.use("Agg")
from pylab import *
+from scipy import stats
import h5py
# Plot parameters
@@ -84,6 +85,24 @@ S = sim["/PartType0/Entropy"][:]
P = sim["/PartType0/Pressure"][:]
rho = sim["/PartType0/Density"][:]
+# Bin te data
+r_bin_edge = np.arange(0., 0.5, 0.01)
+r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
+rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
+v_bin,_,_ = stats.binned_statistic(r, v_r, statistic='mean', bins=r_bin_edge)
+P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
+S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
+u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
+rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
+v2_bin,_,_ = stats.binned_statistic(r, v_r**2, statistic='mean', bins=r_bin_edge)
+P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
+S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
+u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_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)
# Now, work our the solution....
@@ -213,8 +232,9 @@ figure()
# Velocity profile --------------------------------
subplot(231)
-plot(r, v_r, '.', color='r', ms=1.)
+plot(r, v_r, '.', color='r', ms=0.5, alpha=0.2)
plot(r_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Radial~velocity}}~v_r$", labelpad=0)
xlim(0, 1.3 * r_shock)
@@ -222,8 +242,9 @@ ylim(-0.2, 3.8)
# Density profile --------------------------------
subplot(232)
-plot(r, rho, '.', color='r', ms=1.)
+plot(r, rho, '.', color='r', ms=0.5, alpha=0.2)
plot(r_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Density}}~\\rho$", labelpad=2)
xlim(0, 1.3 * r_shock)
@@ -231,8 +252,9 @@ ylim(-0.2, 5.2)
# Pressure profile --------------------------------
subplot(233)
-plot(r, P, '.', color='r', ms=1.)
+plot(r, P, '.', color='r', ms=0.5, alpha=0.2)
plot(r_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Pressure}}~P$", labelpad=0)
xlim(0, 1.3 * r_shock)
@@ -240,8 +262,9 @@ ylim(-1, 12.5)
# Internal energy profile -------------------------
subplot(234)
-plot(r, u, '.', color='r', ms=1.)
+plot(r, u, '.', color='r', ms=0.5, alpha=0.2)
plot(r_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
xlim(0, 1.3 * r_shock)
@@ -249,8 +272,9 @@ ylim(-2, 22)
# Entropy profile ---------------------------------
subplot(235)
-plot(r, S, '.', color='r', ms=1.)
+plot(r, S, '.', color='r', ms=0.5, alpha=0.2)
plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Entropy}}~S$", labelpad=0)
xlim(0, 1.3 * r_shock)
diff --git a/examples/SedovBlast_3D/plotSolution.py b/examples/SedovBlast_3D/plotSolution.py
index 8eca94c98acd7e081cdfc2c785aa1d19d601e81e..6e90a9a43524b3cdb279054764b71fd1b546b366 100644
--- a/examples/SedovBlast_3D/plotSolution.py
+++ b/examples/SedovBlast_3D/plotSolution.py
@@ -35,6 +35,7 @@ gas_gamma = 5./3. # Gas polytropic index
import matplotlib
matplotlib.use("Agg")
from pylab import *
+from scipy import stats
import h5py
# Plot parameters
@@ -85,6 +86,25 @@ S = sim["/PartType0/Entropy"][:]
P = sim["/PartType0/Pressure"][:]
rho = sim["/PartType0/Density"][:]
+# Bin te data
+r_bin_edge = np.arange(0., 0.5, 0.01)
+r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
+rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
+v_bin,_,_ = stats.binned_statistic(r, v_r, statistic='mean', bins=r_bin_edge)
+P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
+S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
+u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
+rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
+v2_bin,_,_ = stats.binned_statistic(r, v_r**2, statistic='mean', bins=r_bin_edge)
+P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
+S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
+u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_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)
+
# Now, work our the solution....
@@ -214,8 +234,9 @@ figure()
# Velocity profile --------------------------------
subplot(231)
-plot(r, v_r, '.', color='r', ms=1.)
+plot(r, v_r, '.', color='r', ms=0.5, alpha=0.2)
plot(r_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Radial~velocity}}~v_r$", labelpad=0)
xlim(0, 1.3 * r_shock)
@@ -223,8 +244,9 @@ ylim(-0.2, 3.8)
# Density profile --------------------------------
subplot(232)
-plot(r, rho, '.', color='r', ms=1.)
+plot(r, rho, '.', color='r', ms=0.5, alpha=0.2)
plot(r_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Density}}~\\rho$", labelpad=2)
xlim(0, 1.3 * r_shock)
@@ -232,8 +254,9 @@ ylim(-0.2, 5.2)
# Pressure profile --------------------------------
subplot(233)
-plot(r, P, '.', color='r', ms=1.)
+plot(r, P, '.', color='r', ms=0.5, alpha=0.2)
plot(r_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Pressure}}~P$", labelpad=0)
xlim(0, 1.3 * r_shock)
@@ -241,8 +264,9 @@ ylim(-1, 12.5)
# Internal energy profile -------------------------
subplot(234)
-plot(r, u, '.', color='r', ms=1.)
+plot(r, u, '.', color='r', ms=0.5, alpha=0.2)
plot(r_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
xlim(0, 1.3 * r_shock)
@@ -250,8 +274,9 @@ ylim(-2, 22)
# Entropy profile ---------------------------------
subplot(235)
-plot(r, S, '.', color='r', ms=1.)
+plot(r, S, '.', color='r', ms=0.5, alpha=0.2)
plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
xlabel("${\\rm{Radius}}~r$", labelpad=0)
ylabel("${\\rm{Entropy}}~S$", labelpad=0)
xlim(0, 1.3 * r_shock)
diff --git a/examples/SodShock_2D/plotSolution.py b/examples/SodShock_2D/plotSolution.py
index 99ba7e9a6e9ae4b6d50688a1428f07e9a08b3b85..b4a203d93518d98ee87282f4ea46d045c4c3b38a 100644
--- a/examples/SodShock_2D/plotSolution.py
+++ b/examples/SodShock_2D/plotSolution.py
@@ -38,6 +38,7 @@ P_R = 0.1 # Pressure right state
import matplotlib
matplotlib.use("Agg")
from pylab import *
+from scipy import stats
import h5py
# Plot parameters
@@ -86,13 +87,30 @@ rho = sim["/PartType0/Density"][:]
N = 1000 # Number of points
x_min = -1.
x_max = 1.
-
x += x_min
-# ---------------------------------------------------------------
-# Don't touch anything after this.
-# ---------------------------------------------------------------
+# 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
@@ -225,8 +243,9 @@ figure()
# Velocity profile --------------------------------
subplot(231)
-plot(x, v, '.', color='r', ms=0.5)
+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)
@@ -234,8 +253,9 @@ ylim(-0.1, 0.95)
# Density profile --------------------------------
subplot(232)
-plot(x, rho, '.', color='r', ms=0.5)
+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)
@@ -243,8 +263,9 @@ ylim(0.05, 1.1)
# Pressure profile --------------------------------
subplot(233)
-plot(x, P, '.', color='r', ms=0.5)
+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)
@@ -252,8 +273,9 @@ ylim(0.01, 1.1)
# Internal energy profile -------------------------
subplot(234)
-plot(x, u, '.', color='r', ms=0.5)
+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)
@@ -261,8 +283,9 @@ ylim(0.8, 2.2)
# Entropy profile ---------------------------------
subplot(235)
-plot(x, S, '.', color='r', ms=0.5)
+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)
diff --git a/examples/SodShock_3D/plotSolution.py b/examples/SodShock_3D/plotSolution.py
index 23a16e6aed73a7281cf78a215940ccdcff722a79..3d9616af55a204db4be9df2e42b355e266944153 100644
--- a/examples/SodShock_3D/plotSolution.py
+++ b/examples/SodShock_3D/plotSolution.py
@@ -38,6 +38,7 @@ P_R = 0.1 # Pressure right state
import matplotlib
matplotlib.use("Agg")
from pylab import *
+from scipy import stats
import h5py
# Plot parameters
@@ -83,16 +84,32 @@ 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
-
-# ---------------------------------------------------------------
-# Don't touch anything after this.
-# ---------------------------------------------------------------
-
+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
@@ -225,8 +242,9 @@ figure()
# Velocity profile --------------------------------
subplot(231)
-plot(x, v, '.', color='r', ms=0.5)
+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)
@@ -234,8 +252,9 @@ ylim(-0.1, 0.95)
# Density profile --------------------------------
subplot(232)
-plot(x, rho, '.', color='r', ms=0.5)
+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)
@@ -243,8 +262,9 @@ ylim(0.05, 1.1)
# Pressure profile --------------------------------
subplot(233)
-plot(x, P, '.', color='r', ms=0.5)
+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)
@@ -252,8 +272,9 @@ ylim(0.01, 1.1)
# Internal energy profile -------------------------
subplot(234)
-plot(x, u, '.', color='r', ms=0.5)
+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)
@@ -261,8 +282,9 @@ ylim(0.8, 2.2)
# Entropy profile ---------------------------------
subplot(235)
-plot(x, S, '.', color='r', ms=0.5)
+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)
diff --git a/examples/Stellar_Disk/makeIC.py b/examples/Stellar_Disk/makeIC.py
deleted file mode 100644
index b5ec65d3fa76d8b377bdaacf3b43a36ab560115e..0000000000000000000000000000000000000000
--- a/examples/Stellar_Disk/makeIC.py
+++ /dev/null
@@ -1,3 +0,0 @@
-import numpy as np
-import h5py as h5
-
diff --git a/examples/Stellar_Disk/stellar_disk.yml b/examples/Stellar_Disk/stellar_disk.yml
deleted file mode 100644
index 1d17f96be7cf476f1959544fd221d7c8b7919915..0000000000000000000000000000000000000000
--- a/examples/Stellar_Disk/stellar_disk.yml
+++ /dev/null
@@ -1 +0,0 @@
-IsothermalPotential
\ No newline at end of file
diff --git a/examples/main.c b/examples/main.c
index 967db3c09eba9e8ae372343e76ddfc3a84c0daf0..852e8d625073bcb97a0367deb93ac5d860768d1d 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -351,7 +351,7 @@ int main(int argc, char *argv[]) {
#endif
/* Initialize unit system and constants */
- struct UnitSystem us;
+ struct unit_system us;
struct phys_const prog_const;
units_init(&us, params, "InternalUnitSystem");
phys_const_init(&us, &prog_const);
diff --git a/src/cell.c b/src/cell.c
index b32ad7c38c892a74d5eddcaa43fbef0a0ef32d58..078310206b50ae65242a3a6b3dac738e2713e8d4 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1271,7 +1271,7 @@ void cell_set_super(struct cell *c, struct cell *super) {
* @param c The #cell.
* @param e The #engine (to get ti_current).
*/
-void cell_drift(struct cell *c, const struct engine *e) {
+void cell_drift_particles(struct cell *c, const struct engine *e) {
const double timeBase = e->timeBase;
const integertime_t ti_old = c->ti_old;
@@ -1295,7 +1295,7 @@ void cell_drift(struct cell *c, const struct engine *e) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
struct cell *cp = c->progeny[k];
- cell_drift(cp, e);
+ cell_drift_particles(cp, e);
dx_max = max(dx_max, cp->dx_max);
h_max = max(h_max, cp->h_max);
}
diff --git a/src/cell.h b/src/cell.h
index aebe773e4ce073dac65cdfc130b2be752e4678d4..c97449f9151b8954ee54cbd6bff5d8ab73cbf53b 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -334,7 +334,7 @@ void cell_check_drift_point(struct cell *c, void *data);
int cell_is_drift_needed(struct cell *c, const struct engine *e);
int cell_unskip_tasks(struct cell *c, struct scheduler *s);
void cell_set_super(struct cell *c, struct cell *super);
-void cell_drift(struct cell *c, const struct engine *e);
+void cell_drift_particles(struct cell *c, const struct engine *e);
void cell_check_timesteps(struct cell *c);
#endif /* SWIFT_CELL_H */
diff --git a/src/common_io.c b/src/common_io.c
index 567ae7249dba788e6b379f8c05eb1139bff6bb94..df0bbdc29ec357da3ba14410c0f9c56e0d69160a 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -134,7 +134,7 @@ int io_is_double_precision(enum IO_DATA_TYPE type) {
*
* @param grp The group from which to read.
* @param name The name of the attribute to read.
- * @param type The #DATA_TYPE of the attribute.
+ * @param type The #IO_DATA_TYPE of the attribute.
* @param data (output) The attribute read from the HDF5 group.
*
* Calls #error() if an error occurs.
@@ -161,7 +161,7 @@ void io_read_attribute(hid_t grp, char* name, enum IO_DATA_TYPE type,
*
* @param grp The group in which to write.
* @param name The name of the attribute to write.
- * @param type The #DATA_TYPE of the attribute.
+ * @param type The #IO_DATA_TYPE of the attribute.
* @param data The attribute to write.
* @param num The number of elements to write
*
@@ -294,11 +294,11 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str) {
/**
* @brief Reads the Unit System from an IC file.
* @param h_file The (opened) HDF5 file from which to read.
- * @param us The UnitSystem to fill.
+ * @param us The unit_system to fill.
*
* If the 'Units' group does not exist in the ICs, cgs units will be assumed
*/
-void io_read_UnitSystem(hid_t h_file, struct UnitSystem* us) {
+void io_read_unit_system(hid_t h_file, struct unit_system* us) {
hid_t h_grp = H5Gopen(h_file, "/Units", H5P_DEFAULT);
@@ -334,11 +334,11 @@ void io_read_UnitSystem(hid_t h_file, struct UnitSystem* us) {
/**
* @brief Writes the current Unit System
* @param h_file The (opened) HDF5 file in which to write
- * @param us The UnitSystem to dump
+ * @param us The unit_system to dump
* @param groupName The name of the HDF5 group to write to
*/
-void io_write_UnitSystem(hid_t h_file, const struct UnitSystem* us,
- const char* groupName) {
+void io_write_unit_system(hid_t h_file, const struct unit_system* us,
+ const char* groupName) {
hid_t h_grpunit = 0;
h_grpunit = H5Gcreate1(h_file, groupName, 0);
diff --git a/src/common_io.h b/src/common_io.h
index 2ebfbbadcfd81ac09275da51a58235fa7e283ee3..577d1664db20196c8d0c48ef5b5e4e960b0e195e 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -68,9 +68,9 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str);
void io_write_code_description(hid_t h_file);
-void io_read_UnitSystem(hid_t h_file, struct UnitSystem* us);
-void io_write_UnitSystem(hid_t h_grp, const struct UnitSystem* us,
- const char* groupName);
+void io_read_unit_system(hid_t h_file, struct unit_system* us);
+void io_write_unit_system(hid_t h_grp, const struct unit_system* us,
+ const char* groupName);
#endif /* defined HDF5 */
diff --git a/src/cooling.c b/src/cooling.c
index e0208dbb591445d0877ef1e703d6e8cf349ddfd6..9452857117367cfc1d4d3eab262c73eba555d4f5 100644
--- a/src/cooling.c
+++ b/src/cooling.c
@@ -34,7 +34,7 @@
* @param cooling The cooling properties to initialize
*/
void cooling_init(const struct swift_params* parameter_file,
- const struct UnitSystem* us,
+ const struct unit_system* us,
const struct phys_const* phys_const,
struct cooling_function_data* cooling) {
diff --git a/src/cooling.h b/src/cooling.h
index f0f227a619182af90d48c9416182f7a2132fd08c..3d50658b9b2cf03cadb6138315b3936d3c4ea4ad 100644
--- a/src/cooling.h
+++ b/src/cooling.h
@@ -42,7 +42,7 @@
/* Common functions */
void cooling_init(const struct swift_params* parameter_file,
- const struct UnitSystem* us,
+ const struct unit_system* us,
const struct phys_const* phys_const,
struct cooling_function_data* cooling);
diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h
index 30ae644bdecbe795794505f64ba1ed767419d82b..fef44e51d78eecbc55a353e809989318b208db50 100644
--- a/src/cooling/const_du/cooling.h
+++ b/src/cooling/const_du/cooling.h
@@ -58,7 +58,7 @@
*/
__attribute__((always_inline)) INLINE static void cooling_cool_part(
const struct phys_const* restrict phys_const,
- const struct UnitSystem* restrict us,
+ const struct unit_system* restrict us,
const struct cooling_function_data* restrict cooling,
struct part* restrict p, struct xpart* restrict xp, float dt) {
@@ -103,7 +103,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
__attribute__((always_inline)) INLINE static float cooling_timestep(
const struct cooling_function_data* restrict cooling,
const struct phys_const* restrict phys_const,
- const struct UnitSystem* restrict us, const struct part* restrict p) {
+ const struct unit_system* restrict us, const struct part* restrict p) {
const float cooling_rate = cooling->cooling_rate;
const float internal_energy = hydro_get_internal_energy(p);
@@ -153,7 +153,7 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
* @param cooling The cooling properties to initialize
*/
static INLINE void cooling_init_backend(
- const struct swift_params* parameter_file, const struct UnitSystem* us,
+ const struct swift_params* parameter_file, const struct unit_system* us,
const struct phys_const* phys_const,
struct cooling_function_data* cooling) {
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index 9fadd51e3c2a3c5462c8476e0aac893e3a2d530d..370d318f255b5c9d6527d288f2e341f37c1e6f40 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -45,7 +45,7 @@
* @param p Pointer to the particle data..
*/
__attribute__((always_inline)) INLINE static float cooling_rate(
- const struct phys_const* const phys_const, const struct UnitSystem* us,
+ const struct phys_const* const phys_const, const struct unit_system* us,
const struct cooling_function_data* cooling, const struct part* p) {
/* Get particle density */
@@ -72,7 +72,7 @@ __attribute__((always_inline)) INLINE static float cooling_rate(
*/
__attribute__((always_inline)) INLINE static void cooling_cool_part(
const struct phys_const* restrict phys_const,
- const struct UnitSystem* restrict us,
+ const struct unit_system* restrict us,
const struct cooling_function_data* restrict cooling,
struct part* restrict p, struct xpart* restrict xp, float dt) {
@@ -112,7 +112,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
__attribute__((always_inline)) INLINE static float cooling_timestep(
const struct cooling_function_data* restrict cooling,
const struct phys_const* restrict phys_const,
- const struct UnitSystem* restrict us, const struct part* restrict p) {
+ const struct unit_system* restrict us, const struct part* restrict p) {
/* Get current internal energy */
const float u = hydro_get_internal_energy(p);
@@ -158,7 +158,7 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
* @param cooling The cooling properties to initialize
*/
static INLINE void cooling_init_backend(
- const struct swift_params* parameter_file, const struct UnitSystem* us,
+ const struct swift_params* parameter_file, const struct unit_system* us,
const struct phys_const* phys_const,
struct cooling_function_data* cooling) {
diff --git a/src/cooling/none/cooling.h b/src/cooling/none/cooling.h
index 0461100dc11e7ffbb4616766923142442b4ac943..63caf27fc6494c6fdd5a52f84eaa371eb824fd16 100644
--- a/src/cooling/none/cooling.h
+++ b/src/cooling/none/cooling.h
@@ -49,7 +49,7 @@
*/
__attribute__((always_inline)) INLINE static void cooling_cool_part(
const struct phys_const* restrict phys_const,
- const struct UnitSystem* restrict us,
+ const struct unit_system* restrict us,
const struct cooling_function_data* restrict cooling,
struct part* restrict p, struct xpart* restrict xp, float dt) {}
@@ -66,7 +66,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
__attribute__((always_inline)) INLINE static float cooling_timestep(
const struct cooling_function_data* restrict cooling,
const struct phys_const* restrict phys_const,
- const struct UnitSystem* restrict us, const struct part* restrict p) {
+ const struct unit_system* restrict us, const struct part* restrict p) {
return FLT_MAX;
}
@@ -107,7 +107,7 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
* @param cooling The cooling properties to initialize
*/
static INLINE void cooling_init_backend(
- const struct swift_params* parameter_file, const struct UnitSystem* us,
+ const struct swift_params* parameter_file, const struct unit_system* us,
const struct phys_const* phys_const,
struct cooling_function_data* cooling) {}
diff --git a/src/engine.c b/src/engine.c
index 12047e81c7c0f20eaa3e6c388a0743a676e3a81d..5f399a4889e8961c59267c7d690f59c6eb6d7a7c 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3365,16 +3365,16 @@ void engine_dump_snapshot(struct engine *e) {
/* Dump... */
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
- write_output_parallel(e, e->snapshotBaseName, e->internalUnits,
+ write_output_parallel(e, e->snapshotBaseName, e->internal_units,
e->snapshotUnits, e->nodeID, e->nr_nodes,
MPI_COMM_WORLD, MPI_INFO_NULL);
#else
- write_output_serial(e, e->snapshotBaseName, e->internalUnits,
+ write_output_serial(e, e->snapshotBaseName, e->internal_units,
e->snapshotUnits, e->nodeID, e->nr_nodes, MPI_COMM_WORLD,
MPI_INFO_NULL);
#endif
#else
- write_output_single(e, e->snapshotBaseName, e->internalUnits,
+ write_output_single(e, e->snapshotBaseName, e->internal_units,
e->snapshotUnits);
#endif
@@ -3466,7 +3466,7 @@ void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
int nr_threads, int with_aff, int policy, int verbose,
enum repartition_type reparttype,
- const struct UnitSystem *internal_units,
+ const struct unit_system *internal_units,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
const struct external_potential *potential,
@@ -3500,7 +3500,7 @@ void engine_init(struct engine *e, struct space *s,
e->timeStep = 0.;
e->timeBase = 0.;
e->timeBase_inv = 0.;
- e->internalUnits = internal_units;
+ e->internal_units = internal_units;
e->timeFirstSnapshot =
parser_get_param_double(params, "Snapshots:time_first");
e->deltaTimeSnapshot =
@@ -3509,7 +3509,7 @@ void engine_init(struct engine *e, struct space *s,
parser_get_param_string(params, "Snapshots:basename", e->snapshotBaseName);
e->snapshotCompression =
parser_get_opt_param_int(params, "Snapshots:compression", 0);
- e->snapshotUnits = malloc(sizeof(struct UnitSystem));
+ e->snapshotUnits = malloc(sizeof(struct unit_system));
units_init_default(e->snapshotUnits, params, "Snapshots", internal_units);
e->dt_min = parser_get_param_double(params, "TimeIntegration:dt_min");
e->dt_max = parser_get_param_double(params, "TimeIntegration:dt_max");
diff --git a/src/engine.h b/src/engine.h
index 50100c687f3a979583d3b3152cf0d7785d0b26d5..633fbaf41de8b2a4afdcf31e87957f4005c921c4 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -141,7 +141,7 @@ struct engine {
size_t updates, g_updates, s_updates;
/* The internal system of units */
- const struct UnitSystem *internalUnits;
+ const struct unit_system *internal_units;
/* Snapshot information */
double timeFirstSnapshot;
@@ -149,7 +149,7 @@ struct engine {
integertime_t ti_nextSnapshot;
char snapshotBaseName[200];
int snapshotCompression;
- struct UnitSystem *snapshotUnits;
+ struct unit_system *snapshotUnits;
/* Statistics information */
FILE *file_stats;
@@ -232,7 +232,7 @@ void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
int nr_threads, int with_aff, int policy, int verbose,
enum repartition_type reparttype,
- const struct UnitSystem *internal_units,
+ const struct unit_system *internal_units,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
const struct external_potential *potential,
diff --git a/src/io_properties.h b/src/io_properties.h
index 4c972ea90c3353c80559a2ee82b9ded571a39d41..9fcf1a1ac67cae6afab6870369e51d06c752fc11 100644
--- a/src/io_properties.h
+++ b/src/io_properties.h
@@ -47,7 +47,7 @@ struct io_props {
enum DATA_IMPORTANCE importance;
/* Units of the quantity */
- enum UnitConversionFactor units;
+ enum unit_conversion_factor units;
/* Pointer to the field of the first particle in the array */
char* field;
@@ -89,7 +89,7 @@ struct io_props {
struct io_props io_make_input_field_(char name[FIELD_BUFFER_SIZE],
enum IO_DATA_TYPE type, int dimension,
enum DATA_IMPORTANCE importance,
- enum UnitConversionFactor units,
+ enum unit_conversion_factor units,
char* field, size_t partSize) {
struct io_props r;
strcpy(r.name, name);
@@ -128,7 +128,7 @@ struct io_props io_make_input_field_(char name[FIELD_BUFFER_SIZE],
*/
struct io_props io_make_output_field_(char name[FIELD_BUFFER_SIZE],
enum IO_DATA_TYPE type, int dimension,
- enum UnitConversionFactor units,
+ enum unit_conversion_factor units,
char* field, size_t partSize) {
struct io_props r;
strcpy(r.name, name);
@@ -171,7 +171,7 @@ struct io_props io_make_output_field_(char name[FIELD_BUFFER_SIZE],
*/
struct io_props io_make_output_field_convert_part_(
char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
- enum UnitConversionFactor units, char* field, size_t partSize,
+ enum unit_conversion_factor units, char* field, size_t partSize,
struct part* parts, float (*functionPtr)(struct engine*, struct part*)) {
struct io_props r;
@@ -215,7 +215,7 @@ struct io_props io_make_output_field_convert_part_(
*/
struct io_props io_make_output_field_convert_gpart_(
char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
- enum UnitConversionFactor units, char* field, size_t partSize,
+ enum unit_conversion_factor units, char* field, size_t partSize,
struct gpart* gparts, float (*functionPtr)(struct engine*, struct gpart*)) {
struct io_props r;
diff --git a/src/parallel_io.c b/src/parallel_io.c
index d210ad05cefd659934e795bfbf9812fef4f82bd7..1d78410fb9163a4c40e6706a91ce6852f7c35a1d 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -71,8 +71,8 @@
*/
void readArray(hid_t grp, const struct io_props prop, size_t N,
long long N_total, long long offset,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* ic_units) {
+ const struct unit_system* internal_units,
+ const struct unit_system* ic_units) {
const size_t typeSize = io_sizeof_type(prop.type);
const size_t copySize = typeSize * prop.dimension;
@@ -191,8 +191,8 @@ void readArray(hid_t grp, const struct io_props prop, size_t N,
* @param N The number of particles to write.
* @param N_total Total number of particles across all cores
* @param offset Offset in the array where this mpi task starts writing
- * @param internal_units The #UnitSystem used internally
- * @param snapshot_units The #UnitSystem used in the snapshots
+ * @param internal_units The #unit_system used internally
+ * @param snapshot_units The #unit_system used in the snapshots
*
* @todo A better version using HDF5 hyper-slabs to write the file directly from
* the part array will be written once the structures have been stabilized.
@@ -201,8 +201,8 @@ void readArray(hid_t grp, const struct io_props prop, size_t N,
void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
char* partTypeGroupName, const struct io_props props, size_t N,
long long N_total, int mpi_rank, long long offset,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* snapshot_units) {
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units) {
const size_t typeSize = io_sizeof_type(props.type);
const size_t copySize = typeSize * props.dimension;
@@ -374,7 +374,7 @@ void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
* Calls #error() if an error occurs.
*
*/
-void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
+void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nstars, int* periodic, int* flag_entropy,
@@ -458,9 +458,9 @@ void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
H5Gclose(h_grp);
/* Read the unit system used in the ICs */
- struct UnitSystem* ic_units = malloc(sizeof(struct UnitSystem));
+ struct unit_system* ic_units = malloc(sizeof(struct unit_system));
if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
- io_read_UnitSystem(h_file, ic_units);
+ io_read_unit_system(h_file, ic_units);
/* Tell the user if a conversion will be needed */
if (mpi_rank == 0) {
@@ -620,8 +620,8 @@ void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
*
* @param e The engine containing all the system.
* @param baseName The common part of the snapshot file name.
- * @param internal_units The #UnitSystem used internally
- * @param snapshot_units The #UnitSystem used in the snapshots
+ * @param internal_units The #unit_system used internally
+ * @param snapshot_units The #unit_system used in the snapshots
* @param mpi_rank The MPI rank of this node.
* @param mpi_size The number of MPI ranks.
* @param comm The MPI communicator.
@@ -636,8 +636,8 @@ void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
*
*/
void write_output_parallel(struct engine* e, const char* baseName,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* snapshot_units,
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units,
int mpi_rank, int mpi_size, MPI_Comm comm,
MPI_Info info) {
@@ -765,10 +765,10 @@ void write_output_parallel(struct engine* e, const char* baseName,
H5Gclose(h_grp);
/* Print the system of Units used in the spashot */
- io_write_UnitSystem(h_file, snapshot_units, "Units");
+ io_write_unit_system(h_file, snapshot_units, "Units");
/* Print the system of Units used internally */
- io_write_UnitSystem(h_file, internal_units, "InternalCodeUnits");
+ io_write_unit_system(h_file, internal_units, "InternalCodeUnits");
/* Tell the user if a conversion will be needed */
if (e->verbose && mpi_rank == 0) {
diff --git a/src/parallel_io.h b/src/parallel_io.h
index e4cb9f5976bc0f5b55207a7422597a05feaa3d5e..d4cde2bd89313f6cf14ce059870925b9ac304b40 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -34,7 +34,7 @@
#if defined(HAVE_HDF5) && defined(WITH_MPI) && defined(HAVE_PARALLEL_HDF5)
-void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
+void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nsparts, int* periodic, int* flag_entropy,
@@ -43,8 +43,8 @@ void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
int dry_run);
void write_output_parallel(struct engine* e, const char* baseName,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* snapshot_units,
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units,
int mpi_rank, int mpi_size, MPI_Comm comm,
MPI_Info info);
#endif
diff --git a/src/physical_constants.c b/src/physical_constants.c
index b726678ec1c0f2c85e409da3898dc3b6b842750f..c851578c96e146261e9512f6899c7b82a8d91097 100644
--- a/src/physical_constants.c
+++ b/src/physical_constants.c
@@ -34,7 +34,8 @@
* @param us The current internal system of units.
* @param internal_const The physical constants to initialize.
*/
-void phys_const_init(struct UnitSystem* us, struct phys_const* internal_const) {
+void phys_const_init(struct unit_system* us,
+ struct phys_const* internal_const) {
/* Units are declared as {U_M, U_L, U_t, U_I, U_T} */
diff --git a/src/physical_constants.h b/src/physical_constants.h
index aac0fa0840e08a1140bbab67e944f7f134f4222b..3731c0ef565c6159592ad2de96a222efc6cf43f2 100644
--- a/src/physical_constants.h
+++ b/src/physical_constants.h
@@ -78,7 +78,7 @@ struct phys_const {
double const_earth_mass;
};
-void phys_const_init(struct UnitSystem* us, struct phys_const* internal_const);
+void phys_const_init(struct unit_system* us, struct phys_const* internal_const);
void phys_const_print(struct phys_const* internal_const);
diff --git a/src/potential.c b/src/potential.c
index 6ee80900952c032d019ad5d20ec086d05d34ef29..94c2a6cc9412d1fc76b70ddebb2edba7878e6209 100644
--- a/src/potential.c
+++ b/src/potential.c
@@ -36,7 +36,7 @@
*/
void potential_init(const struct swift_params* parameter_file,
const struct phys_const* phys_const,
- const struct UnitSystem* us, const struct space* s,
+ const struct unit_system* us, const struct space* s,
struct external_potential* potential) {
potential_init_backend(parameter_file, phys_const, us, s, potential);
diff --git a/src/potential.h b/src/potential.h
index 116ea8302e7f706cdb861540a89d562174d73408..4c20f4c7b4eaa7ca134b3d90ae790a4710a1f217 100644
--- a/src/potential.h
+++ b/src/potential.h
@@ -43,7 +43,7 @@
/* Now, some generic functions, defined in the source file */
void potential_init(const struct swift_params* parameter_file,
const struct phys_const* phys_const,
- const struct UnitSystem* us, const struct space* s,
+ const struct unit_system* us, const struct space* s,
struct external_potential* potential);
void potential_print(const struct external_potential* potential);
diff --git a/src/potential/disc_patch/potential.h b/src/potential/disc_patch/potential.h
index 400539a8d02d29a8d383bb1c523d064f733267c5..8fa40ecd4e6503cde8be00db8c6fb8a70c84ebdf 100644
--- a/src/potential/disc_patch/potential.h
+++ b/src/potential/disc_patch/potential.h
@@ -195,7 +195,7 @@ external_gravity_get_potential_energy(
*/
static INLINE void potential_init_backend(
const struct swift_params* parameter_file,
- const struct phys_const* phys_const, const struct UnitSystem* us,
+ const struct phys_const* phys_const, const struct unit_system* us,
const struct space* s, struct external_potential* potential) {
potential->surface_density = parser_get_param_double(
diff --git a/src/potential/isothermal/potential.h b/src/potential/isothermal/potential.h
index 9c07f3eb67528a003788ca94bd1e2e52dd985a2c..c974618d7b581884f871863bb83200b8cecee7a5 100644
--- a/src/potential/isothermal/potential.h
+++ b/src/potential/isothermal/potential.h
@@ -163,7 +163,7 @@ external_gravity_get_potential_energy(
*/
static INLINE void potential_init_backend(
const struct swift_params* parameter_file,
- const struct phys_const* phys_const, const struct UnitSystem* us,
+ const struct phys_const* phys_const, const struct unit_system* us,
const struct space* s, struct external_potential* potential) {
potential->x =
diff --git a/src/potential/none/potential.h b/src/potential/none/potential.h
index cb6254b4a23b336637cb3c9f36a2dd01170eabad..a8550cad702891ff211539c95c42eca57418c464 100644
--- a/src/potential/none/potential.h
+++ b/src/potential/none/potential.h
@@ -99,7 +99,7 @@ external_gravity_get_potential_energy(
*/
static INLINE void potential_init_backend(
const struct swift_params* parameter_file,
- const struct phys_const* phys_const, const struct UnitSystem* us,
+ const struct phys_const* phys_const, const struct unit_system* us,
const struct space* s, struct external_potential* potential) {}
/**
diff --git a/src/potential/point_mass/potential.h b/src/potential/point_mass/potential.h
index 81b51ab2009ad599d0201708d78c8c64cac991dc..160b74cae82f8dba584e15eac743e6ddc2002c0a 100644
--- a/src/potential/point_mass/potential.h
+++ b/src/potential/point_mass/potential.h
@@ -149,7 +149,7 @@ external_gravity_get_potential_energy(
*/
static INLINE void potential_init_backend(
const struct swift_params* parameter_file,
- const struct phys_const* phys_const, const struct UnitSystem* us,
+ const struct phys_const* phys_const, const struct unit_system* us,
const struct space* s, struct external_potential* potential) {
potential->x =
diff --git a/src/runner.c b/src/runner.c
index b8cb2bcaa37629613a31aab63ca74a45817e7ec7..7cea20de87d8f1602aa1a058cfb8128c8186f680 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -210,7 +210,7 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const struct cooling_function_data *cooling_func = e->cooling_func;
const struct phys_const *constants = e->physical_constants;
- const struct UnitSystem *us = e->internalUnits;
+ const struct unit_system *us = e->internal_units;
const double timeBase = e->timeBase;
TIMER_TIC;
@@ -810,11 +810,11 @@ void runner_do_unskip_mapper(void *map_data, int num_elements,
* @param c The cell.
* @param timer Are we timing this ?
*/
-void runner_do_drift(struct runner *r, struct cell *c, int timer) {
+void runner_do_drift_particles(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
- cell_drift(c, r->e);
+ cell_drift_particles(c, r->e);
if (timer) TIMER_TOC(timer_drift);
}
@@ -834,7 +834,7 @@ void runner_do_drift_mapper(void *map_data, int num_elements,
for (int ind = 0; ind < num_elements; ind++) {
struct cell *c = &cells[ind];
- if (c != NULL && c->nodeID == e->nodeID) cell_drift(c, e);
+ if (c != NULL && c->nodeID == e->nodeID) cell_drift_particles(c, e);
}
}
@@ -1772,7 +1772,7 @@ void *runner_main(void *data) {
break;
#endif
case task_type_drift:
- runner_do_drift(r, ci, 1);
+ runner_do_drift_particles(r, ci, 1);
break;
case task_type_kick1:
runner_do_kick1(r, ci, 1);
diff --git a/src/runner.h b/src/runner.h
index 53e78b00657385c7185e0730d421707c87ccf382..49ad5c98354dac157d30ba292f4e48ce05301586 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -56,7 +56,7 @@ struct runner {
void runner_do_ghost(struct runner *r, struct cell *c, int timer);
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer);
void runner_do_sort(struct runner *r, struct cell *c, int flag, int clock);
-void runner_do_drift(struct runner *r, struct cell *c, int timer);
+void runner_do_drift_particles(struct runner *r, struct cell *c, int timer);
void runner_do_kick1(struct runner *r, struct cell *c, int timer);
void runner_do_kick2(struct runner *r, struct cell *c, int timer);
void runner_do_end_force(struct runner *r, struct cell *c, int timer);
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index 5ebd36ef7bd14564160ff63f356b2e532b09fc22..1fbaa3398dc7cd839ba1b1bafde1a2372658bbf1 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -746,8 +746,8 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
/* Anything to do here? */
if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
- if (!cell_is_drifted(ci, e)) cell_drift(ci, e);
- if (!cell_is_drifted(cj, e)) cell_drift(cj, e);
+ if (!cell_is_drifted(ci, e)) cell_drift_particles(ci, e);
+ if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
/* Get the sort ID. */
double shift[3] = {0.0, 0.0, 0.0};
@@ -1397,7 +1397,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
if (!cell_is_active(c, e)) return;
- if (!cell_is_drifted(c, e)) cell_drift(c, e);
+ if (!cell_is_drifted(c, e)) cell_drift_particles(c, e);
struct part *restrict parts = c->parts;
const int count = c->count;
diff --git a/src/runner_doiact_nosort.h b/src/runner_doiact_nosort.h
index beb7fd7026c23c3a7a516b31d23334095bfddae3..6a442afbdb0c7dfb24c5e9cb386783746a1f05ed 100644
--- a/src/runner_doiact_nosort.h
+++ b/src/runner_doiact_nosort.h
@@ -15,8 +15,8 @@ void DOPAIR1_NOSORT(struct runner *r, struct cell *ci, struct cell *cj) {
/* Anything to do here? */
if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
- if (!cell_is_drifted(ci, e)) cell_drift(ci, e);
- if (!cell_is_drifted(cj, e)) cell_drift(cj, e);
+ if (!cell_is_drifted(ci, e)) cell_drift_particles(ci, e);
+ if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
/* Get the relative distance between the pairs, wrapping. */
double shift[3] = {0.0, 0.0, 0.0};
@@ -140,8 +140,8 @@ void DOPAIR2_NOSORT(struct runner *r, struct cell *ci, struct cell *cj) {
/* Anything to do here? */
if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
- if (!cell_is_drifted(ci, e)) cell_drift(ci, e);
- if (!cell_is_drifted(cj, e)) cell_drift(cj, e);
+ if (!cell_is_drifted(ci, e)) cell_drift_particles(ci, e);
+ if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
/* Get the relative distance between the pairs, wrapping. */
double shift[3] = {0.0, 0.0, 0.0};
diff --git a/src/runner_doiact_vec.c b/src/runner_doiact_vec.c
index b91d288529c0706693d74b0c54d688ee0944aa29..d78de2cd3c508fd17a816316dc931010f2febf80 100644
--- a/src/runner_doiact_vec.c
+++ b/src/runner_doiact_vec.c
@@ -287,7 +287,7 @@ __attribute__((always_inline)) INLINE void runner_doself1_density_vec(
if (!cell_is_active(c, e)) return;
- if (!cell_is_drifted(c, e)) cell_drift(c, e);
+ if (!cell_is_drifted(c, e)) cell_drift_particles(c, e);
/* Get the particle cache from the runner and re-allocate
* the cache if it is not big enough for the cell. */
@@ -536,7 +536,7 @@ __attribute__((always_inline)) INLINE void runner_doself1_density_vec_2(
if (!cell_is_active(c, e)) return;
- if (!cell_is_drifted(c, e)) cell_drift(c, e);
+ if (!cell_is_drifted(c, e)) cell_drift_particles(c, e);
/* TODO: Need to find two active particles, not just one. */
diff --git a/src/serial_io.c b/src/serial_io.c
index 03f7df82c549b49c76799470bab733e0a6820759..28c329070d4995e24a2c375b1967f4e155ea4764 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -74,8 +74,8 @@
*/
void readArray(hid_t grp, const struct io_props props, size_t N,
long long N_total, long long offset,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* ic_units) {
+ const struct unit_system* internal_units,
+ const struct unit_system* ic_units) {
const size_t typeSize = io_sizeof_type(props.type);
const size_t copySize = typeSize * props.dimension;
@@ -182,8 +182,8 @@ void readArray(hid_t grp, const struct io_props props, size_t N,
void prepareArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
char* partTypeGroupName, const struct io_props props,
unsigned long long N_total,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* snapshot_units) {
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units) {
/* Create data space */
const hid_t h_space = H5Screate(H5S_SIMPLE);
@@ -281,14 +281,14 @@ void prepareArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
* @param part_c A (char*) pointer on the first occurrence of the field of
*interest in the parts array
* @param partSize The size in bytes of the particle structure.
- * @param us The UnitSystem currently in use
+ * @param us The unit_system currently in use
* @param convFactor The UnitConversionFactor for this arrayo
*/
void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
char* partTypeGroupName, const struct io_props props, size_t N,
long long N_total, int mpi_rank, long long offset,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* snapshot_units) {
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units) {
const size_t typeSize = io_sizeof_type(props.type);
const size_t copySize = typeSize * props.dimension;
@@ -424,7 +424,7 @@ void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
* @todo Read snapshots distributed in more than one file.
*
*/
-void read_ic_serial(char* fileName, const struct UnitSystem* internal_units,
+void read_ic_serial(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nstars, int* periodic, int* flag_entropy,
@@ -443,7 +443,7 @@ void read_ic_serial(char* fileName, const struct UnitSystem* internal_units,
long long offset[swift_type_count] = {0};
int dimension = 3; /* Assume 3D if nothing is specified */
size_t Ndm = 0;
- struct UnitSystem* ic_units = malloc(sizeof(struct UnitSystem));
+ struct unit_system* ic_units = malloc(sizeof(struct unit_system));
/* First read some information about the content */
if (mpi_rank == 0) {
@@ -507,7 +507,7 @@ void read_ic_serial(char* fileName, const struct UnitSystem* internal_units,
/* Read the unit system used in the ICs */
if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
- io_read_UnitSystem(h_file, ic_units);
+ io_read_unit_system(h_file, ic_units);
if (units_are_equal(ic_units, internal_units)) {
@@ -551,7 +551,7 @@ void read_ic_serial(char* fileName, const struct UnitSystem* internal_units,
MPI_Bcast(periodic, 1, MPI_INT, 0, comm);
MPI_Bcast(&N_total, swift_type_count, MPI_LONG_LONG_INT, 0, comm);
MPI_Bcast(dim, 3, MPI_DOUBLE, 0, comm);
- MPI_Bcast(ic_units, sizeof(struct UnitSystem), MPI_BYTE, 0, comm);
+ MPI_Bcast(ic_units, sizeof(struct unit_system), MPI_BYTE, 0, comm);
/* Divide the particles among the tasks. */
for (int ptype = 0; ptype < swift_type_count; ++ptype) {
@@ -695,8 +695,8 @@ void read_ic_serial(char* fileName, const struct UnitSystem* internal_units,
*
* @param e The engine containing all the system.
* @param baseName The common part of the snapshot file name.
- * @param internal_units The #UnitSystem used internally
- * @param snapshot_units The #UnitSystem used in the snapshots
+ * @param internal_units The #unit_system used internally
+ * @param snapshot_units The #unit_system used in the snapshots
* @param mpi_rank The MPI rank of this node.
* @param mpi_size The number of MPI ranks.
* @param comm The MPI communicator.
@@ -711,8 +711,8 @@ void read_ic_serial(char* fileName, const struct UnitSystem* internal_units,
*
*/
void write_output_serial(struct engine* e, const char* baseName,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* snapshot_units, int mpi_rank,
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units, int mpi_rank,
int mpi_size, MPI_Comm comm, MPI_Info info) {
hid_t h_file = 0, h_grp = 0;
@@ -837,10 +837,10 @@ void write_output_serial(struct engine* e, const char* baseName,
H5Gclose(h_grp);
/* Print the system of Units used in the spashot */
- io_write_UnitSystem(h_file, snapshot_units, "Units");
+ io_write_unit_system(h_file, snapshot_units, "Units");
/* Print the system of Units used internally */
- io_write_UnitSystem(h_file, internal_units, "InternalCodeUnits");
+ io_write_unit_system(h_file, internal_units, "InternalCodeUnits");
/* Tell the user if a conversion will be needed */
if (e->verbose) {
diff --git a/src/serial_io.h b/src/serial_io.h
index 94dd68b93626411ec7dc314d783d80c9e0e967b6..df8d0f1917a69eb28af32aed6780783b0f1099d8 100644
--- a/src/serial_io.h
+++ b/src/serial_io.h
@@ -34,7 +34,7 @@
#if defined(HAVE_HDF5) && defined(WITH_MPI) && !defined(HAVE_PARALLEL_HDF5)
-void read_ic_serial(char* fileName, const struct UnitSystem* internal_units,
+void read_ic_serial(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nstars, int* periodic, int* flag_entropy,
@@ -43,8 +43,8 @@ void read_ic_serial(char* fileName, const struct UnitSystem* internal_units,
int dry_run);
void write_output_serial(struct engine* e, const char* baseName,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* snapshot_units, int mpi_rank,
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units, int mpi_rank,
int mpi_size, MPI_Comm comm, MPI_Info info);
#endif
diff --git a/src/single_io.c b/src/single_io.c
index f8696e96a15078a83069368094cfca198d3bb32a..264e138728458f6a1d19b57e45f4f8f2e6acdc8f 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -59,16 +59,16 @@
* @param h_grp The group from which to read.
* @param prop The #io_props of the field to read
* @param N The number of particles.
- * @param internal_units The #UnitSystem used internally
- * @param ic_units The #UnitSystem used in the ICs
+ * @param internal_units The #unit_system used internally
+ * @param ic_units The #unit_system used in the ICs
*
* @todo A better version using HDF5 hyper-slabs to read the file directly into
*the part array
* will be written once the structures have been stabilized.
*/
void readArray(hid_t h_grp, const struct io_props prop, size_t N,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* ic_units) {
+ const struct unit_system* internal_units,
+ const struct unit_system* ic_units) {
const size_t typeSize = io_sizeof_type(prop.type);
const size_t copySize = typeSize * prop.dimension;
@@ -163,16 +163,16 @@ void readArray(hid_t h_grp, const struct io_props prop, size_t N,
* the HDF5 file.
* @param props The #io_props of the field to read
* @param N The number of particles to write.
- * @param internal_units The #UnitSystem used internally
- * @param snapshot_units The #UnitSystem used in the snapshots
+ * @param internal_units The #unit_system used internally
+ * @param snapshot_units The #unit_system used in the snapshots
*
* @todo A better version using HDF5 hyper-slabs to write the file directly from
* the part array will be written once the structures have been stabilized.
*/
void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
char* partTypeGroupName, const struct io_props props, size_t N,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* snapshot_units) {
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units) {
const size_t typeSize = io_sizeof_type(props.type);
const size_t copySize = typeSize * props.dimension;
@@ -337,7 +337,7 @@ void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
* @todo Read snapshots distributed in more than one file.
*
*/
-void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
+void read_ic_single(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nstars, int* periodic, int* flag_entropy,
@@ -410,9 +410,9 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
H5Gclose(h_grp);
/* Read the unit system used in the ICs */
- struct UnitSystem* ic_units = malloc(sizeof(struct UnitSystem));
+ struct unit_system* ic_units = malloc(sizeof(struct unit_system));
if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
- io_read_UnitSystem(h_file, ic_units);
+ io_read_unit_system(h_file, ic_units);
/* Tell the user if a conversion will be needed */
if (units_are_equal(ic_units, internal_units)) {
@@ -565,8 +565,8 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
*
* @param e The engine containing all the system.
* @param baseName The common part of the snapshot file name.
- * @param internal_units The #UnitSystem used internally
- * @param snapshot_units The #UnitSystem used in the snapshots
+ * @param internal_units The #unit_system used internally
+ * @param snapshot_units The #unit_system used in the snapshots
*
* Creates an HDF5 output file and writes the particles contained
* in the engine. If such a file already exists, it is erased and replaced
@@ -577,8 +577,8 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
*
*/
void write_output_single(struct engine* e, const char* baseName,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* snapshot_units) {
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units) {
hid_t h_file = 0, h_grp = 0;
const size_t Ngas = e->s->nr_parts;
@@ -687,10 +687,10 @@ void write_output_single(struct engine* e, const char* baseName,
H5Gclose(h_grp);
/* Print the system of Units used in the spashot */
- io_write_UnitSystem(h_file, snapshot_units, "Units");
+ io_write_unit_system(h_file, snapshot_units, "Units");
/* Print the system of Units used internally */
- io_write_UnitSystem(h_file, internal_units, "InternalCodeUnits");
+ io_write_unit_system(h_file, internal_units, "InternalCodeUnits");
/* Tell the user if a conversion will be needed */
if (e->verbose) {
diff --git a/src/single_io.h b/src/single_io.h
index bc803b262f70f72ea93090d56112f5a70737c840..e1f1b3a7b39ff3f0d31078cd34b7a6cd5f0dfc77 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -29,7 +29,7 @@
#if defined(HAVE_HDF5) && !defined(WITH_MPI)
-void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
+void read_ic_single(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, size_t* Ngas, size_t* Ndm,
size_t* Nstars, int* periodic, int* flag_entropy,
@@ -37,8 +37,8 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
int dry_run);
void write_output_single(struct engine* e, const char* baseName,
- const struct UnitSystem* internal_units,
- const struct UnitSystem* snapshot_units);
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units);
#endif
diff --git a/src/sourceterms.c b/src/sourceterms.c
index 2a5f326688b52b0e4abfde8dca823d765e0e8a5e..f12071cf912eae3aa8d0e25f0f3b4c5e139de667 100644
--- a/src/sourceterms.c
+++ b/src/sourceterms.c
@@ -37,7 +37,7 @@
* @param source the structure that has all the source term properties
*/
void sourceterms_init(const struct swift_params* parameter_file,
- struct UnitSystem* us, struct sourceterms* source) {
+ struct unit_system* us, struct sourceterms* source) {
#ifdef SOURCETERMS_SN_FEEDBACK
supernova_init(parameter_file, us, source);
#endif /* SOURCETERMS_SN_FEEDBACK */
diff --git a/src/sourceterms.h b/src/sourceterms.h
index 95c37e9bba9881702a6cf6caea1004b8cfb52636..1445bcb777ff634d1e3a2312cb0a49ac155e1020 100644
--- a/src/sourceterms.h
+++ b/src/sourceterms.h
@@ -42,7 +42,7 @@ struct sourceterms {
#endif
void sourceterms_init(const struct swift_params* parameter_file,
- struct UnitSystem* us, struct sourceterms* source);
+ struct unit_system* us, struct sourceterms* source);
void sourceterms_print(struct sourceterms* source);
/**
diff --git a/src/sourceterms/sn_feedback/sn_feedback.h b/src/sourceterms/sn_feedback/sn_feedback.h
index 667c79cb1eef42f7fde79c43059d1baa5c61268e..f2f224ce871ebb768c318aef42a690861dd974df 100644
--- a/src/sourceterms/sn_feedback/sn_feedback.h
+++ b/src/sourceterms/sn_feedback/sn_feedback.h
@@ -171,7 +171,7 @@ __attribute__((always_inline)) INLINE static void supernova_feedback_apply(
*/
__attribute__((always_inline)) INLINE static void supernova_init(
- const struct swift_params* parameter_file, struct UnitSystem* us,
+ const struct swift_params* parameter_file, struct unit_system* us,
struct sourceterms* source) {
source->supernova.time = parser_get_param_double(parameter_file, "SN:time");
source->supernova.energy =
diff --git a/src/timestep.h b/src/timestep.h
index 432f0fd2c4eb713e11272546cfe84e8f6c342cbd..6497e25984ef36d759afd12616b45c8166816dfb 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -109,7 +109,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
float new_dt_cooling = FLT_MAX;
if (e->policy & engine_policy_cooling)
new_dt_cooling = cooling_timestep(e->cooling_func, e->physical_constants,
- e->internalUnits, p);
+ e->internal_units, p);
/* Compute the next timestep (gravity condition) */
float new_dt_grav = FLT_MAX;
diff --git a/src/units.c b/src/units.c
index d50186637f5cc14d8b981d1331c6d656b4575592..13d3e5dd35b14c325527ff35703712258e273ef3 100644
--- a/src/units.c
+++ b/src/units.c
@@ -41,11 +41,11 @@
#include "error.h"
/**
- * @brief Initialises the UnitSystem structure with CGS system
+ * @brief Initialises the unit_system structure with CGS system
*
- * @param us The UnitSystem to initialize
+ * @param us The unit_system to initialize
*/
-void units_init_cgs(struct UnitSystem* us) {
+void units_init_cgs(struct unit_system* us) {
us->UnitMass_in_cgs = 1.;
us->UnitLength_in_cgs = 1.;
@@ -55,14 +55,14 @@ void units_init_cgs(struct UnitSystem* us) {
}
/**
- * @brief Initialises the UnitSystem structure with the constants given in
+ * @brief Initialises the unit_system structure with the constants given in
* rhe parameter file.
*
- * @param us The UnitSystem to initialize.
+ * @param us The unit_system to initialize.
* @param params The parsed parameter file.
* @param category The section of the parameter file to read from.
*/
-void units_init(struct UnitSystem* us, const struct swift_params* params,
+void units_init(struct unit_system* us, const struct swift_params* params,
const char* category) {
char buffer[200];
@@ -80,19 +80,19 @@ void units_init(struct UnitSystem* us, const struct swift_params* params,
}
/**
- * @brief Initialises the UnitSystem structure with the constants given in
+ * @brief Initialises the unit_system structure with the constants given in
* rhe parameter file. Uses a default if the values are not present in the file.
*
- * @param us The UnitSystem to initialize.
+ * @param us The unit_system to initialize.
* @param params The parsed parameter file.
* @param category The section of the parameter file to read from.
* @param def The default unit system to copy from if required.
*/
-void units_init_default(struct UnitSystem* us,
+void units_init_default(struct unit_system* us,
const struct swift_params* params, const char* category,
- const struct UnitSystem* def) {
+ const struct unit_system* def) {
- if (!def) error("Default UnitSystem not allocated");
+ if (!def) error("Default unit_system not allocated");
char buffer[200];
sprintf(buffer, "%s:UnitMass_in_cgs", category);
@@ -116,11 +116,11 @@ void units_init_default(struct UnitSystem* us,
/**
* @brief Returns the base unit conversion factor for a given unit system
- * @param us The UnitSystem used
+ * @param us The unit_system used
* @param baseUnit The base unit
*/
-double units_get_base_unit(const struct UnitSystem* us,
- enum BaseUnits baseUnit) {
+double units_get_base_unit(const struct unit_system* us,
+ enum base_units baseUnit) {
switch (baseUnit) {
case UNIT_MASS:
return us->UnitMass_in_cgs;
@@ -142,7 +142,7 @@ double units_get_base_unit(const struct UnitSystem* us,
* @brief Returns the base unit symbol used internally
* @param baseUnit The base unit
*/
-const char* units_get_base_unit_internal_symbol(enum BaseUnits baseUnit) {
+const char* units_get_base_unit_internal_symbol(enum base_units baseUnit) {
switch (baseUnit) {
case UNIT_MASS:
return "U_M";
@@ -164,7 +164,7 @@ const char* units_get_base_unit_internal_symbol(enum BaseUnits baseUnit) {
* @brief Returns the base unit symbol in the cgs system
* @param baseUnit The base unit
*/
-const char* units_get_base_unit_cgs_symbol(enum BaseUnits baseUnit) {
+const char* units_get_base_unit_cgs_symbol(enum base_units baseUnit) {
switch (baseUnit) {
case UNIT_MASS:
return "g";
@@ -183,7 +183,7 @@ const char* units_get_base_unit_cgs_symbol(enum BaseUnits baseUnit) {
}
void units_get_base_unit_exponants_array(float baseUnitsExp[5],
- enum UnitConversionFactor unit) {
+ enum unit_conversion_factor unit) {
switch (unit) {
case UNIT_CONV_NO_UNITS:
break;
@@ -337,8 +337,8 @@ void units_get_base_unit_exponants_array(float baseUnitsExp[5],
* @param us The system of units in use
* @param unit The unit to convert
*/
-double units_cgs_conversion_factor(const struct UnitSystem* us,
- enum UnitConversionFactor unit) {
+double units_cgs_conversion_factor(const struct unit_system* us,
+ enum unit_conversion_factor unit) {
float baseUnitsExp[5] = {0.f};
units_get_base_unit_exponants_array(baseUnitsExp, unit);
@@ -351,8 +351,8 @@ double units_cgs_conversion_factor(const struct UnitSystem* us,
* @param us The system of units in use
* @param unit The unit to convert
*/
-float units_h_factor(const struct UnitSystem* us,
- enum UnitConversionFactor unit) {
+float units_h_factor(const struct unit_system* us,
+ enum unit_conversion_factor unit) {
float baseUnitsExp[5] = {0.f};
units_get_base_unit_exponants_array(baseUnitsExp, unit);
@@ -365,8 +365,8 @@ float units_h_factor(const struct UnitSystem* us,
* @param us The system of units in use
* @param unit The unit to convert
*/
-float units_a_factor(const struct UnitSystem* us,
- enum UnitConversionFactor unit) {
+float units_a_factor(const struct unit_system* us,
+ enum unit_conversion_factor unit) {
float baseUnitsExp[5] = {0.f};
units_get_base_unit_exponants_array(baseUnitsExp, unit);
@@ -378,8 +378,8 @@ float units_a_factor(const struct UnitSystem* us,
* @brief Returns a string containing the exponents of the base units making up
* the conversion factors
*/
-void units_cgs_conversion_string(char* buffer, const struct UnitSystem* us,
- enum UnitConversionFactor unit) {
+void units_cgs_conversion_string(char* buffer, const struct unit_system* us,
+ enum unit_conversion_factor unit) {
float baseUnitsExp[5] = {0.f};
units_get_base_unit_exponants_array(baseUnitsExp, unit);
@@ -394,13 +394,13 @@ void units_cgs_conversion_string(char* buffer, const struct UnitSystem* us,
* @param baseUnitsExponants The exponent of each base units required to form
* the desired quantity. See conversionFactor() for a working example
*/
-double units_general_cgs_conversion_factor(const struct UnitSystem* us,
+double units_general_cgs_conversion_factor(const struct unit_system* us,
const float baseUnitsExponants[5]) {
double factor = 1.;
for (int i = 0; i < 5; ++i)
if (baseUnitsExponants[i] != 0)
- factor *= pow(units_get_base_unit(us, (enum BaseUnits)i),
+ factor *= pow(units_get_base_unit(us, (enum base_units)i),
baseUnitsExponants[i]);
return factor;
}
@@ -412,7 +412,7 @@ double units_general_cgs_conversion_factor(const struct UnitSystem* us,
* @param baseUnitsExponants The exponent of each base units required to form
* the desired quantity. See conversionFactor() for a working example
*/
-float units_general_h_factor(const struct UnitSystem* us,
+float units_general_h_factor(const struct unit_system* us,
const float baseUnitsExponants[5]) {
float factor_exp = 0.f;
@@ -430,7 +430,7 @@ float units_general_h_factor(const struct UnitSystem* us,
* @param baseUnitsExponants The exponent of each base units required to form
* the desired quantity. See conversionFactor() for a working example
*/
-float units_general_a_factor(const struct UnitSystem* us,
+float units_general_a_factor(const struct unit_system* us,
const float baseUnitsExponants[5]) {
float factor_exp = 0.f;
@@ -449,7 +449,7 @@ float units_general_a_factor(const struct UnitSystem* us,
* the desired quantity. See conversionFactor() for a working example
*/
void units_general_cgs_conversion_string(char* buffer,
- const struct UnitSystem* us,
+ const struct unit_system* us,
const float baseUnitsExponants[5]) {
char temp[14];
const double a_exp = units_general_a_factor(us, baseUnitsExponants);
@@ -493,14 +493,14 @@ void units_general_cgs_conversion_string(char* buffer,
sprintf(temp, " ");
else if (baseUnitsExponants[i] == 1.)
sprintf(temp, "%s ",
- units_get_base_unit_internal_symbol((enum BaseUnits)i));
+ units_get_base_unit_internal_symbol((enum base_units)i));
else if (remainder(baseUnitsExponants[i], 1.) == 0)
sprintf(temp, "%s^%d ",
- units_get_base_unit_internal_symbol((enum BaseUnits)i),
+ units_get_base_unit_internal_symbol((enum base_units)i),
(int)baseUnitsExponants[i]);
else
sprintf(temp, "%s^%7.4f ",
- units_get_base_unit_internal_symbol((enum BaseUnits)i),
+ units_get_base_unit_internal_symbol((enum base_units)i),
baseUnitsExponants[i]);
strncat(buffer, temp, 12);
}
@@ -513,14 +513,15 @@ void units_general_cgs_conversion_string(char* buffer,
if (baseUnitsExponants[i] == 0.)
continue;
else if (baseUnitsExponants[i] == 1.)
- sprintf(temp, "%s ", units_get_base_unit_cgs_symbol((enum BaseUnits)i));
+ sprintf(temp, "%s ",
+ units_get_base_unit_cgs_symbol((enum base_units)i));
else if (remainder(baseUnitsExponants[i], 1.) == 0)
sprintf(temp, "%s^%d ",
- units_get_base_unit_cgs_symbol((enum BaseUnits)i),
+ units_get_base_unit_cgs_symbol((enum base_units)i),
(int)baseUnitsExponants[i]);
else
sprintf(temp, "%s^%7.4f ",
- units_get_base_unit_cgs_symbol((enum BaseUnits)i),
+ units_get_base_unit_cgs_symbol((enum base_units)i),
baseUnitsExponants[i]);
strncat(buffer, temp, 12);
}
@@ -532,11 +533,11 @@ void units_general_cgs_conversion_string(char* buffer,
/**
* @brief Are the two unit systems equal ?
*
- * @param a The First #UnitSystem
- * @param b The second #UnitSystem
+ * @param a The First #unit_system
+ * @param b The second #unit_system
* @return 1 if the systems are the same, 0 otherwise
*/
-int units_are_equal(const struct UnitSystem* a, const struct UnitSystem* b) {
+int units_are_equal(const struct unit_system* a, const struct unit_system* b) {
if (a->UnitMass_in_cgs != b->UnitMass_in_cgs) return 0;
if (a->UnitLength_in_cgs != b->UnitLength_in_cgs) return 0;
@@ -550,13 +551,13 @@ int units_are_equal(const struct UnitSystem* a, const struct UnitSystem* b) {
/**
* @brief Return the unit conversion factor between two systems
*
- * @param from The #UnitSystem we are converting from
- * @param to The #UnitSystem we are converting to
+ * @param from The #unit_system we are converting from
+ * @param to The #unit_system we are converting to
* @param baseUnitsExponants The exponent of each base units required to form
* the desired quantity. See conversionFactor() for a working example
*/
-double units_general_conversion_factor(const struct UnitSystem* from,
- const struct UnitSystem* to,
+double units_general_conversion_factor(const struct unit_system* from,
+ const struct unit_system* to,
const float baseUnitsExponants[5]) {
const double from_cgs =
@@ -570,15 +571,15 @@ double units_general_conversion_factor(const struct UnitSystem* from,
/**
* @brief Return the unit conversion factor between two systems
*
- * @param from The #UnitSystem we are converting from
- * @param to The #UnitSystem we are converting to
+ * @param from The #unit_system we are converting from
+ * @param to The #unit_system we are converting to
* @param unit The unit we are converting
*
* @return The conversion factor
*/
-double units_conversion_factor(const struct UnitSystem* from,
- const struct UnitSystem* to,
- enum UnitConversionFactor unit) {
+double units_conversion_factor(const struct unit_system* from,
+ const struct unit_system* to,
+ enum unit_conversion_factor unit) {
float baseUnitsExp[5] = {0.f};
diff --git a/src/units.h b/src/units.h
index 78fdf1c23c3c276607d5353ee3437d8eb1e96537..a5765495f9f52159ab70a1072c1f8571ddcdf14b 100644
--- a/src/units.h
+++ b/src/units.h
@@ -32,7 +32,7 @@
* internal units. It is used everytime a conversion is performed or an i/o
* function is called.
**/
-struct UnitSystem {
+struct unit_system {
/*! Conversion factor from grams to internal mass units */
double UnitMass_in_cgs;
@@ -54,7 +54,7 @@ struct UnitSystem {
* @brief The base units used in the cgs (and internal) system. All units are
* derived from those.
*/
-enum BaseUnits {
+enum base_units {
UNIT_MASS = 0,
UNIT_LENGTH = 1,
UNIT_TIME = 2,
@@ -65,7 +65,7 @@ enum BaseUnits {
/**
* @brief The different conversion factors supported by default
*/
-enum UnitConversionFactor {
+enum unit_conversion_factor {
UNIT_CONV_NO_UNITS,
UNIT_CONV_MASS,
UNIT_CONV_LENGTH,
@@ -94,47 +94,47 @@ enum UnitConversionFactor {
UNIT_CONV_INV_VOLUME
};
-void units_init_cgs(struct UnitSystem*);
-void units_init(struct UnitSystem*, const struct swift_params*,
+void units_init_cgs(struct unit_system*);
+void units_init(struct unit_system*, const struct swift_params*,
const char* category);
-void units_init_default(struct UnitSystem* us,
+void units_init_default(struct unit_system* us,
const struct swift_params* params, const char* category,
- const struct UnitSystem* def);
+ const struct unit_system* def);
-int units_are_equal(const struct UnitSystem* a, const struct UnitSystem* b);
+int units_are_equal(const struct unit_system* a, const struct unit_system* b);
/* Base units */
-double units_get_base_unit(const struct UnitSystem*, enum BaseUnits);
-const char* units_get_base_unit_internal_symbol(enum BaseUnits);
-const char* units_get_base_unit_cgs_symbol(enum BaseUnits);
+double units_get_base_unit(const struct unit_system*, enum base_units);
+const char* units_get_base_unit_internal_symbol(enum base_units);
+const char* units_get_base_unit_cgs_symbol(enum base_units);
/* Cosmology factors */
-float units_general_h_factor(const struct UnitSystem* us,
+float units_general_h_factor(const struct unit_system* us,
const float baseUnitsExponants[5]);
-float units_h_factor(const struct UnitSystem* us,
- enum UnitConversionFactor unit);
-float units_general_a_factor(const struct UnitSystem* us,
+float units_h_factor(const struct unit_system* us,
+ enum unit_conversion_factor unit);
+float units_general_a_factor(const struct unit_system* us,
const float baseUnitsExponants[5]);
-float units_a_factor(const struct UnitSystem* us,
- enum UnitConversionFactor unit);
+float units_a_factor(const struct unit_system* us,
+ enum unit_conversion_factor unit);
/* Conversion to CGS */
-double units_general_cgs_conversion_factor(const struct UnitSystem* us,
+double units_general_cgs_conversion_factor(const struct unit_system* us,
const float baseUnitsExponants[5]);
-double units_cgs_conversion_factor(const struct UnitSystem* us,
- enum UnitConversionFactor unit);
+double units_cgs_conversion_factor(const struct unit_system* us,
+ enum unit_conversion_factor unit);
void units_general_cgs_conversion_string(char* buffer,
- const struct UnitSystem* us,
+ const struct unit_system* us,
const float baseUnitsExponants[5]);
-void units_cgs_conversion_string(char* buffer, const struct UnitSystem* us,
- enum UnitConversionFactor unit);
+void units_cgs_conversion_string(char* buffer, const struct unit_system* us,
+ enum unit_conversion_factor unit);
/* Conversion between systems */
-double units_general_conversion_factor(const struct UnitSystem* from,
- const struct UnitSystem* to,
+double units_general_conversion_factor(const struct unit_system* from,
+ const struct unit_system* to,
const float baseUnitsExponants[5]);
-double units_conversion_factor(const struct UnitSystem* from,
- const struct UnitSystem* to,
- enum UnitConversionFactor unit);
+double units_conversion_factor(const struct unit_system* from,
+ const struct unit_system* to,
+ enum unit_conversion_factor unit);
#endif /* SWIFT_UNITS_H */
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 91b1cf6dc3b321643aae1f4eec6bd3d7abb48350..aedf0e0b77898cadf97383dfc956a3d6b5ebf39f 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -580,7 +580,7 @@ int main(int argc, char *argv[]) {
// runner_do_kick1(&runner, main_cell, 0);
/* And a gentle drift */
- // runner_do_drift(&runner, main_cell, 0);
+ // runner_do_drift_particles(&runner, main_cell, 0);
/* First, sort stuff */
for (int j = 0; j < 125; ++j) runner_do_sort(&runner, cells[j], 0x1FFF, 0);
@@ -678,7 +678,7 @@ int main(int argc, char *argv[]) {
// runner_do_kick1(&runner, main_cell, 0);
/* And drift it */
- runner_do_drift(&runner, main_cell, 0);
+ runner_do_drift_particles(&runner, main_cell, 0);
/* Initialise the particles */
for (int j = 0; j < 125; ++j) runner_do_init(&runner, cells[j], 0);
diff --git a/tests/testDump.c b/tests/testDump.c
index ab74a1b1f022761efedf5258a20c525fcef47bd6..7343af49f654582de444681ec291311d41251dca 100644
--- a/tests/testDump.c
+++ b/tests/testDump.c
@@ -79,6 +79,9 @@ int main(int argc, char *argv[]) {
/* Finalize the dump. */
dump_close(&d);
+ /* Clean the threads */
+ threadpool_clean(&t);
+
/* Return a happy number. */
return 0;
}
diff --git a/tests/testReading.c b/tests/testReading.c
index cbf25bf880c988bec95a91d5e141bf7554a97fe7..930434e7118ca175185cc2c8f13e83f486e91d2c 100644
--- a/tests/testReading.c
+++ b/tests/testReading.c
@@ -35,7 +35,7 @@ int main() {
struct spart *sparts = NULL;
/* Default unit system */
- struct UnitSystem us;
+ struct unit_system us;
units_init_cgs(&us);
/* Properties of the ICs */