diff --git a/examples/GreshoVortex/makeIC.py b/examples/GreshoVortex/makeIC.py
index 12edcb6e8154ec6f865d28a6daeb02d385d14bbf..b3e0f641293feaf1b6274039fe2bcc82f80b186b 100644
--- a/examples/GreshoVortex/makeIC.py
+++ b/examples/GreshoVortex/makeIC.py
@@ -21,93 +21,83 @@
import h5py
import random
from numpy import *
+import sys
# Generates a swift IC file for the Gresho Vortex in a periodic box
# Parameters
+inputFile = "glass_256.hdf5"
periodic= 1 # 1 For periodic box
-factor = 3
-boxSize = [ 1.0 , 1.0, 1.0/factor ]
-L = 120 # Number of particles along one axis
gamma = 5./3. # Gas adiabatic index
-eta = 1.2349 # 48 ngbs with cubic spline kernel
rho = 1 # Gas density
P0 = 0. # Constant additional pressure (should have no impact on the dynamics)
-fileName = "greshoVortex.hdf5"
-vol = boxSize[0] * boxSize[1] * boxSize[2]
+fileOutputName = "greshoVortex.hdf5"
#---------------------------------------------------
-
-numPart = L**3 / factor
-mass = boxSize[0]*boxSize[1]*boxSize[2] * rho / numPart
-
-#Generate particles
-coords = zeros((numPart, 3))
-v = zeros((numPart, 3))
-m = zeros((numPart, 1))
-h = zeros((numPart, 1))
-u = zeros((numPart, 1))
-ids = zeros((numPart, 1), dtype='L')
-
-partId=0
-for i in range(L):
- for j in range(L):
- for k in range(L/factor):
- index = i*L*L/factor + j*L/factor + k
- x = i * boxSize[0] / L + boxSize[0] / (2*L)
- y = j * boxSize[0] / L + boxSize[0] / (2*L)
- z = k * boxSize[0] / L + boxSize[0] / (2*L)
- r2 = (x - boxSize[0] / 2)**2 + (y - boxSize[1] / 2)**2
- r = sqrt(r2)
- coords[index,0] = x
- coords[index,1] = y
- coords[index,2] = z
- v_phi = 0.
- if r < 0.2:
- v_phi = 5.*r
- elif r < 0.4:
- v_phi = 2. - 5.*r
- else:
- v_phi = 0.
- v[index,0] = -v_phi * (y - boxSize[0] / 2) / r
- v[index,1] = v_phi * (x - boxSize[0] / 2) / r
- v[index,2] = 0.
- m[index] = mass
- h[index] = eta * boxSize[0] / L
- P = P0
- if r < 0.2:
- P = P + 5. + 12.5*r2
- elif r < 0.4:
- P = P + 9. + 12.5*r2 - 20.*r + 4.*log(r/0.2)
- else:
- P = P + 3. + 4.*log(2.)
- u[index] = P / ((gamma - 1.)*rho)
- ids[index] = partId + 1
- partId = partId + 1
-
+fileInput = h5py.File(inputFile, "r")
+header = fileInput["/Header"]
+boxSize = header.attrs["BoxSize"][0]
+coords = fileInput["/PartType0/Coordinates"][:,:]
+h = fileInput["/PartType0/SmoothingLength"][:]
+ids = fileInput["/PartType0/ParticleIDs"][:]
+numPart = size(ids)
+
+m = ones(shape(ids)) * boxSize**2 / numPart
+u = zeros(shape(m))
+v = zeros(shape(coords))
+
+
+for i in range(numPart):
+ x = coords[i,0]
+ y = coords[i,1]
+ z = coords[i,2]
+
+ r2 = (x - boxSize / 2)**2 + (y - boxSize / 2)**2
+ r = sqrt(r2)
+
+ v_phi = 0.
+ if r < 0.2:
+ v_phi = 5.*r
+ elif r < 0.4:
+ v_phi = 2. - 5.*r
+ else:
+ v_phi = 0.
+ v[i,0] = -v_phi * (y - boxSize / 2) / r
+ v[i,1] = v_phi * (x - boxSize / 2) / r
+ v[i,2] = 0.
+
+ P = P0
+ if r < 0.2:
+ P = P + 5. + 12.5*r2
+ elif r < 0.4:
+ P = P + 9. + 12.5*r2 - 20.*r + 4.*log(r/0.2)
+ else:
+ P = P + 3. + 4.*log(2.)
+ u[i] = P / ((gamma - 1.)*rho)
+
#File
-file = h5py.File(fileName, 'w')
+fileOutput = h5py.File(fileOutputName, 'w')
# Header
-grp = file.create_group("/Header")
+grp = fileOutput.create_group("/Header")
grp.attrs["BoxSize"] = boxSize
grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0
-grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["NumFileOutputsPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
#Runtime parameters
-grp = file.create_group("/RuntimePars")
+grp = fileOutput.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
#Units
-grp = file.create_group("/Units")
+grp = fileOutput.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.
@@ -115,20 +105,20 @@ 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 = fileOutput.create_group("/PartType0")
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
ds[()] = coords
ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
-ds = grp.create_dataset('Masses', (numPart,1), 'f')
-ds[()] = m
+ds = grp.create_dataset('Masses', (numPart, 1), 'f')
+ds[()] = m.reshape((numPart,1))
ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-ds[()] = h
+ds[()] = h.reshape((numPart,1))
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-ds[()] = u
+ds[()] = u.reshape((numPart,1))
ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
-ds[()] = ids[:]
+ds[()] = ids.reshape((numPart,1))
-file.close()
+fileOutput.close()
diff --git a/examples/PerturbedPlane/makeIC.py b/examples/PerturbedPlane/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..20b720419ff095016daad23828b81ca880ea9c2e
--- /dev/null
+++ b/examples/PerturbedPlane/makeIC.py
@@ -0,0 +1,115 @@
+###############################################################################
+ # 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
+import sys
+import random
+from numpy import *
+
+# Generates a swift IC file containing a perturbed cartesian distribution of particles
+# at a constant density and pressure in a cubic box
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 1.
+L = int(sys.argv[1]) # Number of particles along one axis
+rho = 1. # Density
+P = 1. # Pressure
+gamma = 5./3. # Gas adiabatic index
+pert = 0.1 # Perturbation scale (in units of the interparticle separation)
+fileName = "perturbedPlane.hdf5"
+
+
+#---------------------------------------------------
+numPart = L**2
+mass = boxSize**2 * rho / numPart
+internalEnergy = P / ((gamma - 1.)*rho)
+
+#Generate particles
+coords = zeros((numPart, 3))
+v = zeros((numPart, 3))
+m = zeros((numPart, 1))
+h = zeros((numPart, 1))
+u = zeros((numPart, 1))
+ids = zeros((numPart, 1), dtype='L')
+
+for i in range(L):
+ for j in range(L):
+ index = i*L + j
+ x = i * boxSize / L + boxSize / (2*L) + random.random() * pert * boxSize/(2.*L)
+ y = j * boxSize / L + boxSize / (2*L) + random.random() * pert * boxSize/(2.*L)
+ z = 0
+ coords[index,0] = x
+ coords[index,1] = y
+ coords[index,2] = z
+ v[index,0] = 0.
+ v[index,1] = 0.
+ v[index,2] = 0.
+ m[index] = mass
+ h[index] = 1.23485 * boxSize / L
+ u[index] = internalEnergy
+ ids[index] = index
+
+
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total"] = numPart
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+#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")
+ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = coords
+ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+ds = grp.create_dataset('Masses', (numPart,1), 'f')
+ds[()] = m
+ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
+ds[()] = h
+ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
+ds[()] = u
+ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
+ds[()] = ids + 1
+
+file.close()
diff --git a/examples/PerturbedPlane/perturbedPlane.yml b/examples/PerturbedPlane/perturbedPlane.yml
new file mode 100644
index 0000000000000000000000000000000000000000..e810131cc4b1c66b46b483b1605f9d84bcf203b3
--- /dev/null
+++ b/examples/PerturbedPlane/perturbedPlane.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: 10. # The end time of the simulation (in internal units).
+ dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: perturbedPlane # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1e-1 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-3 # 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.
+ max_smoothing_length: 0.1 # Maximal smoothing length allowed (in internal units).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./perturbedPlane.hdf5 # The file to read
diff --git a/src/const.h b/src/const.h
index 6710fe9e52dfdbede468282d84d4441b923a59fe..d9e4e09d2081361778df329660db3b6ad8f2e204 100644
--- a/src/const.h
+++ b/src/const.h
@@ -37,8 +37,8 @@
#define const_max_u_change 0.1f
/* Dimensionality of the problem */
-#define HYDRO_DIMENSION_3D
-//#define HYDRO_DIMENSION_2D
+//#define HYDRO_DIMENSION_3D
+#define HYDRO_DIMENSION_2D
//#define HYDRO_DIMENSION_1D
/* Hydrodynamical adiabatic index. */