diff --git a/.gitignore b/.gitignore
index 3521b353ec4b1c4dea7192047e197596f083fe99..1a43373acd789366119becb30662be2855db4a51 100644
--- a/.gitignore
+++ b/.gitignore
@@ -29,6 +29,7 @@ examples/used_parameters.yml
examples/energy.txt
examples/*/*.xmf
examples/*/*.hdf5
+examples/*/*.png
examples/*/*.txt
examples/*/used_parameters.yml
examples/*/*/*.xmf
diff --git a/examples/GreshoVortex/getGlass.sh b/examples/GreshoVortex/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ae3c977064f5e7a408aa249c5fd9089b3c52ecb1
--- /dev/null
+++ b/examples/GreshoVortex/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/GreshoVortex/gresho.yml b/examples/GreshoVortex/gresho.yml
new file mode 100644
index 0000000000000000000000000000000000000000..066b433a948db232c9665fd7d007753a19e058fa
--- /dev/null
+++ b/examples/GreshoVortex/gresho.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: 1. # 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: gresho # 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-2 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ 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: ./greshoVortex.hdf5 # The file to read
diff --git a/examples/GreshoVortex/makeIC.py b/examples/GreshoVortex/makeIC.py
index 12edcb6e8154ec6f865d28a6daeb02d385d14bbf..2f5bebc00ce0f86d3f4f3cccd030cfff5f90d51d 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
+# Generates a swift IC file for the Gresho-Chan vortex in a periodic box
# Parameters
-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
+rho0 = 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"
+fileGlass = "glassPlane_128.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
-
+# Get position and smoothing lengths from the glass
+fileInput = h5py.File(fileGlass, 'r')
+coords = fileInput["/PartType0/Coordinates"][:,:]
+h = fileInput["/PartType0/SmoothingLength"][:]
+ids = fileInput["/PartType0/ParticleIDs"][:]
+boxSize = fileInput["/Header"].attrs["BoxSize"][0]
+numPart = size(h)
+fileInput.close()
+
+# Now generate the rest
+m = ones(numPart) * rho0 * boxSize**2 / numPart
+u = zeros(numPart)
+v = zeros((numPart, 3))
+
+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.)*rho0)
+
#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.attrs["PeriodicBoundariesOn"] = periodic
+grp = fileOutput.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
#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/GreshoVortex/plotSolution.py b/examples/GreshoVortex/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..e9d957d80cff3db61d7095d63b46c122c98526ce
--- /dev/null
+++ b/examples/GreshoVortex/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 Gresho-Chan vortex and plots the SPH answer
+
+# Parameters
+gas_gamma = 5./3. # Gas adiabatic index
+rho0 = 1 # Gas density
+P0 = 0. # Constant additional pressure (should have no impact on the dynamics)
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+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])
+
+# Generate the analytic solution at this time
+N = 200
+R_max = 0.8
+solution_r = arange(0, R_max, R_max / N)
+solution_P = zeros(N)
+solution_v_phi = zeros(N)
+solution_v_r = zeros(N)
+
+for i in range(N):
+ if solution_r[i] < 0.2:
+ solution_P[i] = P0 + 5. + 12.5*solution_r[i]**2
+ solution_v_phi[i] = 5.*solution_r[i]
+ elif solution_r[i] < 0.4:
+ solution_P[i] = P0 + 9. + 12.5*solution_r[i]**2 - 20.*solution_r[i] + 4.*log(solution_r[i]/0.2)
+ solution_v_phi[i] = 2. -5.*solution_r[i]
+ else:
+ solution_P[i] = P0 + 3. + 4.*log(2.)
+ solution_v_phi[i] = 0.
+
+solution_rho = ones(N) * rho0
+solution_s = solution_P / solution_rho**gas_gamma
+solution_u = solution_P /((gas_gamma - 1.)*solution_rho)
+
+# Read the simulation data
+sim = h5py.File("gresho_%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"]
+
+pos = sim["/PartType0/Coordinates"][:,:]
+x = pos[:,0] - boxSize / 2
+y = pos[:, 1] - boxSize / 2
+vel = sim["/PartType0/Velocities"][:,:]
+r = sqrt(x**2 + y**2)
+v_r = (x * vel[:,0] + y * vel[:,1]) / r
+v_phi = (-y * vel[:,0] + x * vel[:,1]) / r
+v_norm = sqrt(vel[:,0]**2 + vel[:,1]**2)
+rho = sim["/PartType0/Density"][:]
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+
+# Plot the interesting quantities
+figure()
+
+
+# Azimuthal velocity profile -----------------------------
+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)
+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)
+ylabel("${\\rm{Azimuthal~velocity}}~v_\\phi$", labelpad=0)
+xlim(0,R_max)
+ylim(-0.1, 1.2)
+
+# Radial density profile --------------------------------
+subplot(232)
+
+plot(r, rho, '.', color='r', ms=0.5)
+plot(solution_r, solution_rho, '--', color='k', alpha=0.8, 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)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+xlim(0,R_max)
+ylim(rho0-0.3, rho0 + 0.3)
+#yticks([-0.2, -0.1, 0., 0.1, 0.2])
+
+# Radial pressure profile --------------------------------
+subplot(233)
+
+plot(r, P, '.', color='r', ms=0.5)
+plot(solution_r, solution_P, '--', color='k', alpha=0.8, 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)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(0, R_max)
+ylim(4.9 + P0, P0 + 6.1)
+
+# Internal energy profile --------------------------------
+subplot(234)
+
+plot(r, u, '.', color='r', ms=0.5)
+plot(solution_r, solution_u, '--', color='k', alpha=0.8, 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)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(0,R_max)
+ylim(7.3, 9.1)
+
+
+# Radial entropy profile --------------------------------
+subplot(235)
+
+plot(r, S, '.', color='r', ms=0.5)
+plot(solution_r, solution_s, '--', color='k', alpha=0.8, 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)
+ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+xlim(0, R_max)
+ylim(4.9 + P0, P0 + 6.1)
+
+# Image --------------------------------------------------
+#subplot(234)
+#scatter(pos[:,0], pos[:,1], c=v_norm, cmap="PuBu", edgecolors='face', s=4, vmin=0, vmax=1)
+#text(0.95, 0.95, "$|v|$", ha="right", va="top")
+#xlim(0,1)
+#ylim(0,1)
+#xlabel("$x$", labelpad=0)
+#ylabel("$y$", labelpad=0)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Gresho-Chan vortex with $\\gamma=%.3f$ at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "Background $\\rho_0=%.3f$"%rho0, fontsize=10)
+text(-0.49, 0.7, "Background $P_0=%.3f$"%P0, 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("GreshoVortex.png", dpi=200)
diff --git a/examples/GreshoVortex/run.sh b/examples/GreshoVortex/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..69e28350a2cb9e6747eb879288603f968ebf0a40
--- /dev/null
+++ b/examples/GreshoVortex/run.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e glass_128.hdf5 ]
+then
+ echo "Fetching initial glass file for the Gresho-Chan vortex example..."
+ ./getGlass.sh
+fi
+if [ ! -e greshoVortex.hdf5 ]
+then
+ echo "Generating initial conditions for the Gresho-Chan vortex example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 1 gresho.yml
+
+# Plot the solution
+python plotSolution.py 11
diff --git a/examples/GreshoVortex/solution.py b/examples/GreshoVortex/solution.py
deleted file mode 100644
index 5b282caf2d7f3311ddb595781fed74c51bb4819f..0000000000000000000000000000000000000000
--- a/examples/GreshoVortex/solution.py
+++ /dev/null
@@ -1,69 +0,0 @@
-###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- # 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 random
-from numpy import *
-
-# Computes the analytical solution of the Gresho-Chan vortex
-# The script works for a given initial box and background pressure and computes the solution for any time t (The solution is constant over time).
-# The code writes five files rho.dat, P.dat, v.dat, u.dat and s.dat with the density, pressure, internal energy and
-# entropic function on N radial points between r=0 and r=R_max.
-
-# Parameters
-rho0 = 1. # Background Density
-P0 = 0. # Background Pressure
-gamma = 5./3. # Gas polytropic index
-N = 1000 # Number of radial points
-R_max = 1. # Maximal radius
-
-# ---------------------------------------------------------------
-# Don't touch anything after this.
-# ---------------------------------------------------------------
-
-r = arange(0, R_max, R_max / N)
-rho = ones(N)
-P = zeros(N)
-v = zeros(N)
-u = zeros(N)
-s = zeros(N)
-
-
-for i in range(N):
- if r[i] < 0.2:
- P[i] = P0 + 5. + 12.5*r[i]**2
- v[i] = 5.*r[i]
- elif r[i] < 0.4:
- P[i] = P0 + 9. + 12.5*r[i]**2 - 20.*r[i] + 4.*log(r[i]/0.2)
- v[i] = 2. -5.*r[i]
- else:
- P[i] = P0 + 3. + 4.*log(2.)
- v[i] = 0.
- rho[i] = rho0
- s[i] = P[i] / rho[i]**gamma
- u[i] = P[i] /((gamma - 1.)*rho[i])
-
-savetxt("rho.dat", column_stack((r, rho)))
-savetxt("P.dat", column_stack((r, P)))
-savetxt("v.dat", column_stack((r, v)))
-savetxt("u.dat", column_stack((r, u)))
-savetxt("s.dat", column_stack((r, s)))
-
-
-
diff --git a/examples/PerturbedBox_2D/makeIC.py b/examples/PerturbedBox_2D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..20b720419ff095016daad23828b81ca880ea9c2e
--- /dev/null
+++ b/examples/PerturbedBox_2D/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/PerturbedBox_2D/perturbedPlane.yml b/examples/PerturbedBox_2D/perturbedPlane.yml
new file mode 100644
index 0000000000000000000000000000000000000000..e810131cc4b1c66b46b483b1605f9d84bcf203b3
--- /dev/null
+++ b/examples/PerturbedBox_2D/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/examples/PerturbedBox/makeIC.py b/examples/PerturbedBox_3D/makeIC.py
similarity index 100%
rename from examples/PerturbedBox/makeIC.py
rename to examples/PerturbedBox_3D/makeIC.py
diff --git a/examples/PerturbedBox/perturbedBox.yml b/examples/PerturbedBox_3D/perturbedBox.yml
similarity index 100%
rename from examples/PerturbedBox/perturbedBox.yml
rename to examples/PerturbedBox_3D/perturbedBox.yml
diff --git a/examples/PerturbedBox/run.sh b/examples/PerturbedBox_3D/run.sh
similarity index 100%
rename from examples/PerturbedBox/run.sh
rename to examples/PerturbedBox_3D/run.sh
diff --git a/examples/SedovBlast/makeIC.py b/examples/SedovBlast_3D/makeIC.py
similarity index 100%
rename from examples/SedovBlast/makeIC.py
rename to examples/SedovBlast_3D/makeIC.py
diff --git a/examples/SedovBlast/makeIC_fcc.py b/examples/SedovBlast_3D/makeIC_fcc.py
similarity index 100%
rename from examples/SedovBlast/makeIC_fcc.py
rename to examples/SedovBlast_3D/makeIC_fcc.py
diff --git a/examples/SedovBlast/profile.py b/examples/SedovBlast_3D/profile.py
similarity index 100%
rename from examples/SedovBlast/profile.py
rename to examples/SedovBlast_3D/profile.py
diff --git a/examples/SedovBlast/rdf.py b/examples/SedovBlast_3D/rdf.py
similarity index 100%
rename from examples/SedovBlast/rdf.py
rename to examples/SedovBlast_3D/rdf.py
diff --git a/examples/SedovBlast/run.sh b/examples/SedovBlast_3D/run.sh
similarity index 100%
rename from examples/SedovBlast/run.sh
rename to examples/SedovBlast_3D/run.sh
diff --git a/examples/SedovBlast/sedov.py b/examples/SedovBlast_3D/sedov.py
similarity index 100%
rename from examples/SedovBlast/sedov.py
rename to examples/SedovBlast_3D/sedov.py
diff --git a/examples/SedovBlast/sedov.yml b/examples/SedovBlast_3D/sedov.yml
similarity index 100%
rename from examples/SedovBlast/sedov.yml
rename to examples/SedovBlast_3D/sedov.yml
diff --git a/examples/SedovBlast/solution.py b/examples/SedovBlast_3D/solution.py
similarity index 100%
rename from examples/SedovBlast/solution.py
rename to examples/SedovBlast_3D/solution.py
diff --git a/examples/SodShock_1D/makeIC.py b/examples/SodShock_1D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..e024188e867c1f2636187c2c53157c214752d6f7
--- /dev/null
+++ b/examples/SodShock_1D/makeIC.py
@@ -0,0 +1,118 @@
+###############################################################################
+ # 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 Sod Shock in a periodic box
+
+# Parameters
+gamma = 5./3. # Gas adiabatic index
+numPart_L = 800 # Number of particles in the left state
+x_min = -1.
+x_max = 1.
+rho_L = 1. # Density left state
+rho_R = 0.125 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
+fileName = "sodShock.hdf5"
+
+
+#---------------------------------------------------
+
+# Find how many particles we actually have
+boxSize = x_max - x_min
+numPart_R = int(numPart_L * (rho_R / rho_L))
+numPart = numPart_L + numPart_R
+
+# Now get the distances
+delta_L = (boxSize/2) / numPart_L
+delta_R = (boxSize/2) / numPart_R
+offset_L = delta_L / 2
+offset_R = delta_R / 2
+
+# Build the arrays
+coords = zeros((numPart, 3))
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+h = zeros(numPart)
+u = zeros(numPart)
+
+# Set the particles on the left
+for i in range(numPart_L):
+ coords[i,0] = x_min + offset_L + i * delta_L
+ u[i] = P_L / (rho_L * (gamma - 1.))
+ h[i] = 1.2348 * delta_L
+ m[i] = boxSize * rho_L / (2. * numPart_L)
+ v[i,0] = v_L
+
+# Set the particles on the right
+for j in range(numPart_R):
+ i = numPart_L + j
+ coords[i,0] = offset_R + j * delta_R
+ u[i] = P_R / (rho_R * (gamma - 1.))
+ h[i] = 1.2348 * delta_R
+ m[i] = boxSize * rho_R / (2. * numPart_R)
+ v[i,0] = v_R
+
+# Shift particles
+coords[:,0] -= x_min
+
+#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
+
+#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/SodShock_1D/plotSolution.py b/examples/SodShock_1D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..57f66fe29b9be86a3c4de0d90eafe615d0cb2dbb
--- /dev/null
+++ b/examples/SodShock_1D/plotSolution.py
@@ -0,0 +1,290 @@
+###############################################################################
+ # 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 Sod shock and plots the SPH answer
+
+
+# Generates the analytical solution for the Sod shock test case
+# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
+# The code writes five files rho.dat, P.dat, v.dat, u.dat and s.dat with the density, pressure, internal energy and
+# entropic function on N points between x_min and x_max.
+# This follows the solution given in (Toro, 2009)
+
+
+# Parameters
+gas_gamma = 5./3. # Polytropic index
+rho_L = 1. # Density left state
+rho_R = 0.125 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
+
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left' : 0.045,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.05,
+'figure.subplot.top' : 0.99,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+snap = int(sys.argv[1])
+
+
+# Read the simulation data
+sim = h5py.File("sodShock_%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 = 1000 # Number of points
+x_min = -1.
+x_max = 1.
+
+x += x_min
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+c_L = sqrt(gas_gamma * P_L / rho_L) # Speed of the rarefaction wave
+c_R = sqrt(gas_gamma * P_R / rho_R) # Speed of the shock front
+
+# Helpful variable
+Gama = (gas_gamma - 1.) / (gas_gamma + 1.)
+beta = (gas_gamma - 1.) / (2. * gas_gamma)
+
+# Characteristic function and its derivative, following Toro (2009)
+def compute_f(P_3, P, c):
+ u = P_3 / P
+ if u > 1:
+ term1 = gas_gamma*((gas_gamma+1.)*u + gas_gamma-1.)
+ term2 = sqrt(2./term1)
+ fp = (u - 1.)*c*term2
+ dfdp = c*term2/P + (u - 1.)*c/term2*(-1./term1**2)*gas_gamma*(gas_gamma+1.)/P
+ else:
+ fp = (u**beta - 1.)*(2.*c/(gas_gamma-1.))
+ dfdp = 2.*c/(gas_gamma-1.)*beta*u**(beta-1.)/P
+ return (fp, dfdp)
+
+# Solution of the Riemann problem following Toro (2009)
+def RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R):
+ P_new = ((c_L + c_R + (v_L - v_R)*0.5*(gas_gamma-1.))/(c_L / P_L**beta + c_R / P_R**beta))**(1./beta)
+ P_3 = 0.5*(P_R + P_L)
+ f_L = 1.
+ while fabs(P_3 - P_new) > 1e-6:
+ P_3 = P_new
+ (f_L, dfdp_L) = compute_f(P_3, P_L, c_L)
+ (f_R, dfdp_R) = compute_f(P_3, P_R, c_R)
+ f = f_L + f_R + (v_R - v_L)
+ df = dfdp_L + dfdp_R
+ dp = -f/df
+ prnew = P_3 + dp
+ v_3 = v_L - f_L
+ return (P_new, v_3)
+
+
+# Solve Riemann problem for post-shock region
+(P_3, v_3) = RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R)
+
+# Check direction of shocks and wave
+shock_R = (P_3 > P_R)
+shock_L = (P_3 > P_L)
+
+# Velocity of shock front and and rarefaction wave
+if shock_R:
+ v_right = v_R + c_R**2*(P_3/P_R - 1.)/(gas_gamma*(v_3-v_R))
+else:
+ v_right = c_R + 0.5*(gas_gamma+1.)*v_3 - 0.5*(gas_gamma-1.)*v_R
+
+if shock_L:
+ v_left = v_L + c_L**2*(P_3/p_L - 1.)/(gas_gamma*(v_3-v_L))
+else:
+ v_left = c_L - 0.5*(gas_gamma+1.)*v_3 + 0.5*(gas_gamma-1.)*v_L
+
+# Compute position of the transitions
+x_23 = -fabs(v_left) * time
+if shock_L :
+ x_12 = -fabs(v_left) * time
+else:
+ x_12 = -(c_L - v_L) * time
+
+x_34 = v_3 * time
+
+x_45 = fabs(v_right) * time
+if shock_R:
+ x_56 = fabs(v_right) * time
+else:
+ x_56 = (c_R + v_R) * time
+
+
+# Prepare arrays
+delta_x = (x_max - x_min) / N
+x_s = arange(x_min, x_max, delta_x)
+rho_s = zeros(N)
+P_s = zeros(N)
+v_s = zeros(N)
+
+# Compute solution in the different regions
+for i in range(N):
+ if x_s[i] <= x_12:
+ rho_s[i] = rho_L
+ P_s[i] = P_L
+ v_s[i] = v_L
+ if x_s[i] >= x_12 and x_s[i] < x_23:
+ if shock_L:
+ rho_s[i] = rho_L*(Gama + P_3/P_L)/(1. + Gama * P_3/P_L)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ else:
+ rho_s[i] = rho_L*(Gama * (0. - x_s[i])/(c_L * time) + Gama * v_L/c_L + (1.-Gama))**(2./(gas_gamma-1.))
+ P_s[i] = P_L*(rho_s[i] / rho_L)**gas_gamma
+ v_s[i] = (1.-Gama)*(c_L -(0. - x_s[i]) / time) + Gama*v_L
+ if x_s[i] >= x_23 and x_s[i] < x_34:
+ if shock_L:
+ rho_s[i] = rho_L*(Gama + P_3/P_L)/(1+Gama * P_3/p_L)
+ else:
+ rho_s[i] = rho_L*(P_3 / P_L)**(1./gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_34 and x_s[i] < x_45:
+ if shock_R:
+ rho_s[i] = rho_R*(Gama + P_3/P_R)/(1. + Gama * P_3/P_R)
+ else:
+ rho_s[i] = rho_R*(P_3 / P_R)**(1./gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_45 and x_s[i] < x_56:
+ if shock_R:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+ else:
+ rho_s[i] = rho_R*(Gama*(x_s[i])/(c_R*time) - Gama*v_R/c_R + (1.-Gama))**(2./(gas_gamma-1.))
+ P_s[i] = p_R*(rho_s[i]/rho_R)**gas_gamma
+ v_s[i] = (1.-Gama)*(-c_R - (-x_s[i])/time) + Gama*v_R
+ if x_s[i] >= x_56:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+
+
+# Additional arrays
+u_s = P_s / (rho_s * (gas_gamma - 1.)) #internal energy
+s_s = P_s / rho_s**gas_gamma # entropic function
+
+
+# Plot the interesting quantities
+figure()
+
+# Velocity profile --------------------------------
+subplot(231)
+plot(x, v, '.', color='r', ms=4.0)
+plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(-0.1, 0.95)
+
+# Density profile --------------------------------
+subplot(232)
+plot(x, rho, '.', color='r', ms=4.0)
+plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.05, 1.1)
+
+# Pressure profile --------------------------------
+subplot(233)
+plot(x, P, '.', color='r', ms=4.0)
+plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.01, 1.1)
+
+# Internal energy profile -------------------------
+subplot(234)
+plot(x, u, '.', color='r', ms=4.0)
+plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 2.2)
+
+# Entropy 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=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 3.8)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Sod shock with $\\gamma=%.3f$ in 1D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "Left:~~ $(P_L, \\rho_L, v_L) = (%.3f, %.3f, %.3f)$"%(P_L, rho_L, v_L), fontsize=10)
+text(-0.49, 0.7, "Right: $(P_R, \\rho_R, v_R) = (%.3f, %.3f, %.3f)$"%(P_R, rho_R, v_R), fontsize=10)
+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("SodShock.png", dpi=200)
diff --git a/examples/SodShock_1D/run.sh b/examples/SodShock_1D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..e3ac218c56caa81d0e7a6817e03a8db20bb575d5
--- /dev/null
+++ b/examples/SodShock_1D/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e sodShock.hdf5 ]
+then
+ echo "Generating initial conditions for the 1D SodShock example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 1 sodShock.yml
+
+# Plot the result
+python plotSolution.py 1
diff --git a/examples/SodShock_1D/sodShock.yml b/examples/SodShock_1D/sodShock.yml
new file mode 100644
index 0000000000000000000000000000000000000000..d5c4d0b034ff5351222d2162e37e3e40ceab834f
--- /dev/null
+++ b/examples/SodShock_1D/sodShock.yml
@@ -0,0 +1,36 @@
+# 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.2 # 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-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sodShock # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.2 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ max_smoothing_length: 0.4 # 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: ./sodShock.hdf5 # The file to read
+
diff --git a/examples/SodShock_2D/getGlass.sh b/examples/SodShock_2D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..f4cb4ebcb4b452b2b123462bc97eed532f43ba25
--- /dev/null
+++ b/examples/SodShock_2D/getGlass.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_128.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_48.hdf5
diff --git a/examples/SodShock_2D/makeIC.py b/examples/SodShock_2D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..ac2b9ab45fb68921bce7971c46048b344955140d
--- /dev/null
+++ b/examples/SodShock_2D/makeIC.py
@@ -0,0 +1,122 @@
+###############################################################################
+ # 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 Sod Shock in a periodic box
+
+# Parameters
+gamma = 5./3. # Gas adiabatic index
+x_min = -1.
+x_max = 1.
+rho_L = 1. # Density left state
+rho_R = 0.140625 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
+fileName = "sodShock.hdf5"
+
+
+#---------------------------------------------------
+boxSize = (x_max - x_min)
+
+glass_L = h5py.File("glassPlane_128.hdf5", "r")
+glass_R = h5py.File("glassPlane_48.hdf5", "r")
+
+pos_L = glass_L["/PartType0/Coordinates"][:,:] * 0.5
+pos_R = glass_R["/PartType0/Coordinates"][:,:] * 0.5
+h_L = glass_L["/PartType0/SmoothingLength"][:] * 0.5
+h_R = glass_R["/PartType0/SmoothingLength"][:] * 0.5
+
+# Merge things
+aa = pos_L - array([0.5, 0., 0.])
+pos_LL = append(pos_L, pos_L + array([0.5, 0., 0.]), axis=0)
+pos_RR = append(pos_R, pos_R + array([0.5, 0., 0.]), axis=0)
+pos = append(pos_LL - array([1.0, 0., 0.]), pos_RR, axis=0)
+h_LL = append(h_L, h_L)
+h_RR = append(h_R, h_R)
+h = append(h_LL, h_RR)
+
+numPart_L = size(h_LL)
+numPart_R = size(h_RR)
+numPart = size(h)
+
+vol_L = 0.5
+vol_R = 0.5
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+u = zeros(numPart)
+
+for i in range(numPart):
+ x = pos[i,0]
+
+ if x < 0: #left
+ u[i] = P_L / (rho_L * (gamma - 1.))
+ m[i] = rho_L * vol_L / numPart_L
+ v[i,0] = v_L
+ else: #right
+ u[i] = P_R / (rho_R * (gamma - 1.))
+ m[i] = rho_R * vol_R / numPart_R
+ v[i,0] = v_R
+
+# Shift particles
+pos[:,0] -= x_min
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [boxSize, 0.5, 0.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
+
+#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/SodShock_2D/plotSolution.py b/examples/SodShock_2D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..e74651c567bb711b3190662cf78d421a66134775
--- /dev/null
+++ b/examples/SodShock_2D/plotSolution.py
@@ -0,0 +1,290 @@
+###############################################################################
+ # 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 Sod shock and plots the SPH answer
+
+
+# Generates the analytical solution for the Sod shock test case
+# The script works for a given left (x<0) and right (x>0) state and computes the solution at a later time t.
+# The code writes five files rho.dat, P.dat, v.dat, u.dat and s.dat with the density, pressure, internal energy and
+# entropic function on N points between x_min and x_max.
+# This follows the solution given in (Toro, 2009)
+
+
+# Parameters
+gas_gamma = 5./3. # Polytropic index
+rho_L = 1. # Density left state
+rho_R = 0.140625 # Density right state
+v_L = 0. # Velocity left state
+v_R = 0. # Velocity right state
+P_L = 1. # Pressure left state
+P_R = 0.1 # Pressure right state
+
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left' : 0.045,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.05,
+'figure.subplot.top' : 0.99,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+snap = int(sys.argv[1])
+
+
+# Read the simulation data
+sim = h5py.File("sodShock_%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 = 1000 # Number of points
+x_min = -1.
+x_max = 1.
+
+x += x_min
+
+# ---------------------------------------------------------------
+# Don't touch anything after this.
+# ---------------------------------------------------------------
+
+c_L = sqrt(gas_gamma * P_L / rho_L) # Speed of the rarefaction wave
+c_R = sqrt(gas_gamma * P_R / rho_R) # Speed of the shock front
+
+# Helpful variable
+Gama = (gas_gamma - 1.) / (gas_gamma + 1.)
+beta = (gas_gamma - 1.) / (2. * gas_gamma)
+
+# Characteristic function and its derivative, following Toro (2009)
+def compute_f(P_3, P, c):
+ u = P_3 / P
+ if u > 1:
+ term1 = gas_gamma*((gas_gamma+1.)*u + gas_gamma-1.)
+ term2 = sqrt(2./term1)
+ fp = (u - 1.)*c*term2
+ dfdp = c*term2/P + (u - 1.)*c/term2*(-1./term1**2)*gas_gamma*(gas_gamma+1.)/P
+ else:
+ fp = (u**beta - 1.)*(2.*c/(gas_gamma-1.))
+ dfdp = 2.*c/(gas_gamma-1.)*beta*u**(beta-1.)/P
+ return (fp, dfdp)
+
+# Solution of the Riemann problem following Toro (2009)
+def RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R):
+ P_new = ((c_L + c_R + (v_L - v_R)*0.5*(gas_gamma-1.))/(c_L / P_L**beta + c_R / P_R**beta))**(1./beta)
+ P_3 = 0.5*(P_R + P_L)
+ f_L = 1.
+ while fabs(P_3 - P_new) > 1e-6:
+ P_3 = P_new
+ (f_L, dfdp_L) = compute_f(P_3, P_L, c_L)
+ (f_R, dfdp_R) = compute_f(P_3, P_R, c_R)
+ f = f_L + f_R + (v_R - v_L)
+ df = dfdp_L + dfdp_R
+ dp = -f/df
+ prnew = P_3 + dp
+ v_3 = v_L - f_L
+ return (P_new, v_3)
+
+
+# Solve Riemann problem for post-shock region
+(P_3, v_3) = RiemannProblem(rho_L, P_L, v_L, rho_R, P_R, v_R)
+
+# Check direction of shocks and wave
+shock_R = (P_3 > P_R)
+shock_L = (P_3 > P_L)
+
+# Velocity of shock front and and rarefaction wave
+if shock_R:
+ v_right = v_R + c_R**2*(P_3/P_R - 1.)/(gas_gamma*(v_3-v_R))
+else:
+ v_right = c_R + 0.5*(gas_gamma+1.)*v_3 - 0.5*(gas_gamma-1.)*v_R
+
+if shock_L:
+ v_left = v_L + c_L**2*(P_3/p_L - 1.)/(gas_gamma*(v_3-v_L))
+else:
+ v_left = c_L - 0.5*(gas_gamma+1.)*v_3 + 0.5*(gas_gamma-1.)*v_L
+
+# Compute position of the transitions
+x_23 = -fabs(v_left) * time
+if shock_L :
+ x_12 = -fabs(v_left) * time
+else:
+ x_12 = -(c_L - v_L) * time
+
+x_34 = v_3 * time
+
+x_45 = fabs(v_right) * time
+if shock_R:
+ x_56 = fabs(v_right) * time
+else:
+ x_56 = (c_R + v_R) * time
+
+
+# Prepare arrays
+delta_x = (x_max - x_min) / N
+x_s = arange(x_min, x_max, delta_x)
+rho_s = zeros(N)
+P_s = zeros(N)
+v_s = zeros(N)
+
+# Compute solution in the different regions
+for i in range(N):
+ if x_s[i] <= x_12:
+ rho_s[i] = rho_L
+ P_s[i] = P_L
+ v_s[i] = v_L
+ if x_s[i] >= x_12 and x_s[i] < x_23:
+ if shock_L:
+ rho_s[i] = rho_L*(Gama + P_3/P_L)/(1. + Gama * P_3/P_L)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ else:
+ rho_s[i] = rho_L*(Gama * (0. - x_s[i])/(c_L * time) + Gama * v_L/c_L + (1.-Gama))**(2./(gas_gamma-1.))
+ P_s[i] = P_L*(rho_s[i] / rho_L)**gas_gamma
+ v_s[i] = (1.-Gama)*(c_L -(0. - x_s[i]) / time) + Gama*v_L
+ if x_s[i] >= x_23 and x_s[i] < x_34:
+ if shock_L:
+ rho_s[i] = rho_L*(Gama + P_3/P_L)/(1+Gama * P_3/p_L)
+ else:
+ rho_s[i] = rho_L*(P_3 / P_L)**(1./gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_34 and x_s[i] < x_45:
+ if shock_R:
+ rho_s[i] = rho_R*(Gama + P_3/P_R)/(1. + Gama * P_3/P_R)
+ else:
+ rho_s[i] = rho_R*(P_3 / P_R)**(1./gas_gamma)
+ P_s[i] = P_3
+ v_s[i] = v_3
+ if x_s[i] >= x_45 and x_s[i] < x_56:
+ if shock_R:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+ else:
+ rho_s[i] = rho_R*(Gama*(x_s[i])/(c_R*time) - Gama*v_R/c_R + (1.-Gama))**(2./(gas_gamma-1.))
+ P_s[i] = p_R*(rho_s[i]/rho_R)**gas_gamma
+ v_s[i] = (1.-Gama)*(-c_R - (-x_s[i])/time) + Gama*v_R
+ if x_s[i] >= x_56:
+ rho_s[i] = rho_R
+ P_s[i] = P_R
+ v_s[i] = v_R
+
+
+# Additional arrays
+u_s = P_s / (rho_s * (gas_gamma - 1.)) #internal energy
+s_s = P_s / rho_s**gas_gamma # entropic function
+
+
+# Plot the interesting quantities
+figure()
+
+# Velocity profile --------------------------------
+subplot(231)
+plot(x, v, '.', color='r', ms=0.5)
+plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(-0.1, 0.95)
+
+# Density profile --------------------------------
+subplot(232)
+plot(x, rho, '.', color='r', ms=0.5)
+plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.05, 1.1)
+
+# Pressure profile --------------------------------
+subplot(233)
+plot(x, P, '.', color='r', ms=0.5)
+plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.01, 1.1)
+
+# Internal energy profile -------------------------
+subplot(234)
+plot(x, u, '.', color='r', ms=0.5)
+plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 2.2)
+
+# Entropy profile ---------------------------------
+subplot(235)
+plot(x, S, '.', color='r', ms=0.5)
+plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+xlim(-0.5, 0.5)
+ylim(0.8, 3.8)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Sod shock with $\\gamma=%.3f$ in 2D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "Left:~~ $(P_L, \\rho_L, v_L) = (%.3f, %.3f, %.3f)$"%(P_L, rho_L, v_L), fontsize=10)
+text(-0.49, 0.7, "Right: $(P_R, \\rho_R, v_R) = (%.3f, %.3f, %.3f)$"%(P_R, rho_R, v_R), fontsize=10)
+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("SodShock.png", dpi=200)
diff --git a/examples/SodShock_2D/run.sh b/examples/SodShock_2D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..3d1f4beb3ac4c4becfda24ed39390efec8074da7
--- /dev/null
+++ b/examples/SodShock_2D/run.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e glassPlane_128.hdf5 ]
+then
+ echo "Fetching initial glass file for the Sod shock example..."
+ ./getGlass.sh
+fi
+if [ ! -e sodShock.hdf5 ]
+then
+ echo "Generating initial conditions for the Sod shock example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 1 sodShock.yml
+
+python plotSolution.py 1
diff --git a/examples/SodShock_2D/sodShock.yml b/examples/SodShock_2D/sodShock.yml
new file mode 100644
index 0000000000000000000000000000000000000000..6805724ff58defebc41f3fb5b636d0003b0d6680
--- /dev/null
+++ b/examples/SodShock_2D/sodShock.yml
@@ -0,0 +1,36 @@
+# 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.2 # 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-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sodShock # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.2 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ max_smoothing_length: 0.02 # 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: ./sodShock.hdf5 # The file to read
+
diff --git a/examples/SodShock/glass_001.hdf5 b/examples/SodShock_3D/glass_001.hdf5
similarity index 100%
rename from examples/SodShock/glass_001.hdf5
rename to examples/SodShock_3D/glass_001.hdf5
diff --git a/examples/SodShock/glass_002.hdf5 b/examples/SodShock_3D/glass_002.hdf5
similarity index 100%
rename from examples/SodShock/glass_002.hdf5
rename to examples/SodShock_3D/glass_002.hdf5
diff --git a/examples/SodShock/makeIC.py b/examples/SodShock_3D/makeIC.py
similarity index 100%
rename from examples/SodShock/makeIC.py
rename to examples/SodShock_3D/makeIC.py
diff --git a/examples/SodShock/rhox.py b/examples/SodShock_3D/rhox.py
similarity index 100%
rename from examples/SodShock/rhox.py
rename to examples/SodShock_3D/rhox.py
diff --git a/examples/SodShock/run.sh b/examples/SodShock_3D/run.sh
similarity index 100%
rename from examples/SodShock/run.sh
rename to examples/SodShock_3D/run.sh
diff --git a/examples/SodShock/sodShock.yml b/examples/SodShock_3D/sodShock.yml
similarity index 100%
rename from examples/SodShock/sodShock.yml
rename to examples/SodShock_3D/sodShock.yml
diff --git a/examples/SodShock/solution.py b/examples/SodShock_3D/solution.py
similarity index 100%
rename from examples/SodShock/solution.py
rename to examples/SodShock_3D/solution.py
diff --git a/examples/UniformBox_2D/makeIC.py b/examples/UniformBox_2D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..41b7d695a2376b990703706977ef111be8f3a355
--- /dev/null
+++ b/examples/UniformBox_2D/makeIC.py
@@ -0,0 +1,114 @@
+###############################################################################
+ # 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
+fileName = "uniformPlane.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) + boxSize/(2.*L)
+ y = j * boxSize / L + boxSize / (2*L) + 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/UniformBox_2D/run.sh b/examples/UniformBox_2D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..fbabedb1cd564e0831228591ce32d9a73fe08d9a
--- /dev/null
+++ b/examples/UniformBox_2D/run.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e uniformPlane.hdf5 ]
+then
+ echo "Generating initial conditions for the uniform box example..."
+ python makeIC.py 100
+fi
+
+../swift -s -t 16 uniformPlane.yml
diff --git a/examples/UniformBox_2D/uniformPlane.yml b/examples/UniformBox_2D/uniformPlane.yml
new file mode 100644
index 0000000000000000000000000000000000000000..a3e2d275e50fb20f66ea6494c1202319e462dbed
--- /dev/null
+++ b/examples/UniformBox_2D/uniformPlane.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: uniformPlane # 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: ./uniformPlane.hdf5 # The file to read
diff --git a/examples/UniformBox/makeIC.py b/examples/UniformBox_3D/makeIC.py
similarity index 100%
rename from examples/UniformBox/makeIC.py
rename to examples/UniformBox_3D/makeIC.py
diff --git a/examples/UniformBox/makeICbig.py b/examples/UniformBox_3D/makeICbig.py
similarity index 100%
rename from examples/UniformBox/makeICbig.py
rename to examples/UniformBox_3D/makeICbig.py
diff --git a/examples/UniformBox/run.sh b/examples/UniformBox_3D/run.sh
similarity index 100%
rename from examples/UniformBox/run.sh
rename to examples/UniformBox_3D/run.sh
diff --git a/examples/UniformBox/uniformBox.yml b/examples/UniformBox_3D/uniformBox.yml
similarity index 100%
rename from examples/UniformBox/uniformBox.yml
rename to examples/UniformBox_3D/uniformBox.yml
diff --git a/src/Makefile.am b/src/Makefile.am
index 7e2ceeb6067bb831620472b50ccd51fdfd50b494..4d6a67c7569464486a017d8306c3e54730ebb3b2 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -56,7 +56,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \
units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
- equation_of_state.h \
+ dimension.h equation_of_state.h \
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h
index 34d7095d6c6bc7c82794b60ebbcf6dfc48db033f..78d1cc9d2deaf60c5d933b74a089c77446b44414 100644
--- a/src/adiabatic_index.h
+++ b/src/adiabatic_index.h
@@ -19,6 +19,12 @@
#ifndef SWIFT_ADIABATIC_INDEX_H
#define SWIFT_ADIABATIC_INDEX_H
+/**
+ * @file adiabatic_index.h
+ * @brief Defines the adiabatic index (polytropix index) \f$\gamma\f$ of the
+ * problem and (fast) mathematical functions involving it.
+ */
+
/* Config parameters. */
#include "../config.h"
diff --git a/src/const.h b/src/const.h
index 498b60b711e88e8054b5e0e73a18475d98277448..6710fe9e52dfdbede468282d84d4441b923a59fe 100644
--- a/src/const.h
+++ b/src/const.h
@@ -36,6 +36,11 @@
/* Time integration constants. */
#define const_max_u_change 0.1f
+/* Dimensionality of the problem */
+#define HYDRO_DIMENSION_3D
+//#define HYDRO_DIMENSION_2D
+//#define HYDRO_DIMENSION_1D
+
/* Hydrodynamical adiabatic index. */
#define HYDRO_GAMMA_5_3
//#define HYDRO_GAMMA_4_3
diff --git a/src/dimension.h b/src/dimension.h
new file mode 100644
index 0000000000000000000000000000000000000000..6395d4d04e50d40b733e7a74dafb7d0ab277d204
--- /dev/null
+++ b/src/dimension.h
@@ -0,0 +1,180 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DIMENSION_H
+#define SWIFT_DIMENSION_H
+
+/**
+ * @file dimension.h
+ * @brief Defines the dimensionality \f$d\f$ of the problem and (fast)
+ * mathematical functions involving it
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "const.h"
+#include "inline.h"
+#include "vector.h"
+
+/* First define some constants */
+#if defined(HYDRO_DIMENSION_3D)
+
+#define hydro_dimension 3.f
+#define hydro_dimension_inv 0.3333333333f
+#define hydro_dimention_unit_sphere ((float)(4. * M_PI / 3.))
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+#define hydro_dimension 2.f
+#define hydro_dimension_inv 0.5f
+#define hydro_dimention_unit_sphere ((float)M_PI)
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+#define hydro_dimension 1.f
+#define hydro_dimension_inv 1.f
+#define hydro_dimention_unit_sphere 2.f
+
+#else
+
+#error "A problem dimensionality must be chosen in const.h !"
+
+#endif
+
+/**
+ * @brief Returns the argument to the power given by the dimension
+ *
+ * Computes \f$x^d\f$.
+ */
+__attribute__((always_inline)) INLINE static float pow_dimension(float x) {
+
+#if defined(HYDRO_DIMENSION_3D)
+
+ return x * x * x;
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+ return x * x;
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+ return x;
+
+#else
+
+ error("The dimension is not defined !");
+ return 0.f;
+
+#endif
+}
+
+/**
+ * @brief Returns the argument to the power given by the dimension plus one
+ *
+ * Computes \f$x^{d+1}\f$.
+ */
+__attribute__((always_inline)) INLINE static float pow_dimension_plus_one(
+ float x) {
+
+#if defined(HYDRO_DIMENSION_3D)
+
+ const float x2 = x * x;
+ return x2 * x2;
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+ return x * x * x;
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+ return x * x;
+
+#else
+
+ error("The dimension is not defined !");
+ return 0.f;
+
+#endif
+}
+
+/* ------------------------------------------------------------------------- */
+#ifdef WITH_VECTORIZATION
+
+/**
+ * @brief Returns the argument to the power given by the dimension (vector
+ * version)
+ *
+ * Computes \f$x^d\f$.
+ */
+__attribute__((always_inline)) INLINE static vector pow_dimension_vec(
+ vector x) {
+
+#if defined(HYDRO_DIMENSION_3D)
+
+ return (vector)(x.v * x.v * x.v);
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+ return (vector)(x.v * x.v);
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+ return x;
+
+#else
+
+ error("The dimension is not defined !");
+ return vec_set(0.f);
+
+#endif
+}
+
+/**
+ * @brief Returns the argument to the power given by the dimension plus one
+ * (vector version)
+ *
+ * Computes \f$x^{d+1}\f$.
+ */
+__attribute__((always_inline)) INLINE static vector pow_dimension_plus_one_vec(
+ vector x) {
+
+#if defined(HYDRO_DIMENSION_3D)
+
+ const vector x2 = (vector)(x.v * x.v);
+ return (vector)(x2.v * x2.v);
+
+#elif defined(HYDRO_DIMENSION_2D)
+
+ return (vector)(x.v * x.v * x.v);
+
+#elif defined(HYDRO_DIMENSION_1D)
+
+ return (vector)(x.v * x.v);
+
+#else
+
+ error("The dimension is not defined !");
+ return vec_set(0.f);
+
+#endif
+}
+#endif
+
+#endif /* SWIFT_DIMENSION_H */
diff --git a/src/drift.h b/src/drift.h
index a6e347cc337385f89bece53bc477bff7d163c202..880595dc59e3e5174ee5e888e595a9204ad383e2 100644
--- a/src/drift.h
+++ b/src/drift.h
@@ -25,6 +25,7 @@
/* Local headers. */
#include "const.h"
#include "debug.h"
+#include "dimension.h"
#include "hydro.h"
#include "part.h"
@@ -85,7 +86,7 @@ __attribute__((always_inline)) INLINE static void drift_part(
p->h *= expf(w1);
/* Predict density */
- const float w2 = -3.0f * w1;
+ const float w2 = -hydro_dimension * w1;
if (fabsf(w2) < 0.2f)
p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
else
diff --git a/src/engine.c b/src/engine.c
index ebd589bbb0af01cec3196c658c4d014daca8df38..a980bcbc9c2f42127fb43b30c7d5573a8e9957ef 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2411,17 +2411,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
struct clocks_time time1, time2;
clocks_gettime(&time1);
- if (e->nodeID == 0) message("Initialising particles");
-
- /* Make sure all particles are ready to go */
- /* i.e. clean-up any stupid state in the ICs */
- if (e->policy & engine_policy_hydro) {
- space_map_cells_pre(s, 0, cell_init_parts, NULL);
- }
- if ((e->policy & engine_policy_self_gravity) ||
- (e->policy & engine_policy_external_gravity)) {
- space_map_cells_pre(s, 0, cell_init_gparts, NULL);
- }
+ if (e->nodeID == 0) message("Running initialisation fake time-step.");
engine_prepare(e);
diff --git a/src/equation_of_state.h b/src/equation_of_state.h
index d9f1a46f240a8bdf965749a5513ca3448d08aeb9..b4a36e8a3ef0bcc281d1f939f89fde08ecf00be9 100644
--- a/src/equation_of_state.h
+++ b/src/equation_of_state.h
@@ -19,6 +19,13 @@
#ifndef SWIFT_EQUATION_OF_STATE_H
#define SWIFT_EQUATION_OF_STATE_H
+/**
+ * @file equation_of_state.h
+ * @brief Defines the equation of state of the gas we simulate in the form of
+ * relations between thermodynamic quantities. These are later used internally
+ * by all hydro schemes
+ */
+
/* Config parameters. */
#include "../config.h"
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 5903ce1384f616a123c74b579539fe53747d17e4..755d4be527e4a81b2b4d2b3b829ea16a74ccc7c5 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -81,6 +81,9 @@ gravity_compute_timestep_self(const struct phys_const* const phys_const,
*/
__attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
struct gpart* gp) {
+
+ gp->ti_begin = 0;
+ gp->ti_end = 0;
gp->epsilon = 0.; // MATTHIEU
}
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 7fe78f75fc00569ee8d34ea30a822bac394531e2..51e09a5d943a7f3782e484289faddece2567f15d 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -137,7 +137,15 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
- struct part *restrict p, struct xpart *restrict xp) {}
+ struct part *restrict p, struct xpart *restrict xp) {
+
+ p->ti_begin = 0;
+ p->ti_end = 0;
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
+ xp->u_full = p->u;
+}
/**
* @brief Prepares a particle for the density calculation.
@@ -173,33 +181,33 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Some smoothing length multiples. */
const float h = p->h;
- const float ih = 1.0f / h;
- const float ih2 = ih * ih;
- const float ih4 = ih2 * ih2;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
/* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root;
- p->rho_dh -= 3.0f * p->mass * kernel_root;
+ p->rho_dh -= hydro_dimension * p->mass * kernel_root;
p->density.wcount += kernel_root;
/* Finish the calculation by inserting the missing h-factors */
- p->rho *= ih * ih2;
- p->rho_dh *= ih4;
+ p->rho *= h_inv_dim;
+ p->rho_dh *= h_inv_dim_plus_one;
p->density.wcount *= kernel_norm;
- p->density.wcount_dh *= ih * kernel_gamma * kernel_norm;
+ p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
const float irho = 1.f / p->rho;
/* Compute the derivative term */
- p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
+ p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
/* Finish calculation of the velocity curl components */
- p->density.rot_v[0] *= ih4 * irho;
- p->density.rot_v[1] *= ih4 * irho;
- p->density.rot_v[2] *= ih4 * irho;
+ p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
+ p->density.rot_v[1] *= h_inv_dim_plus_one * irho;
+ p->density.rot_v[2] *= h_inv_dim_plus_one * irho;
/* Finish calculation of the velocity divergence */
- p->density.div_v *= ih4 * irho;
+ p->density.div_v *= h_inv_dim_plus_one * irho;
}
/**
@@ -217,7 +225,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Some smoothing length multiples. */
const float h = p->h;
- const float ih = 1.0f / h;
+ const float h_inv = 1.0f / h;
/* Pre-compute some stuff for the balsara switch. */
const float normDiv_v = fabs(p->density.div_v);
@@ -231,11 +239,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
sqrtf(hydro_gamma * hydro_gamma_minus_one * u);
/* Compute the P/Omega/rho2. */
- xp->omega = 1.0f + 0.3333333333f * h * p->rho_dh / p->rho;
+ xp->omega = 1.0f + hydro_dimension_inv * h * p->rho_dh / p->rho;
p->force.P_over_rho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega);
/* Balsara switch */
- p->force.balsara = normDiv_v / (normDiv_v + normRot_v + 0.0001f * fc * ih);
+ p->force.balsara = normDiv_v / (normDiv_v + normRot_v + 0.0001f * fc * h_inv);
/* Viscosity parameter decay time */
const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
@@ -309,7 +317,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part *restrict p) {
- p->force.h_dt *= p->h * 0.333333333f;
+ p->force.h_dt *= p->h * hydro_dimension_inv;
}
/**
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index 141ab8f884da7f120f2b636933904f0d795a5d8b..4975cac35aaa3d4b4b9f99319dfd061998014453 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -75,7 +75,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
kernel_deval(xi, &wi, &wi_dx);
pi->rho += mj * wi;
- pi->rho_dh -= mj * (3.0f * wi + xi * wi_dx);
+ pi->rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= xi * wi_dx;
@@ -88,7 +88,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
kernel_deval(xj, &wj, &wj_dx);
pj->rho += mi * wj;
- pj->rho_dh -= mi * (3.0f * wj + xj * wj_dx);
+ pj->rho_dh -= mi * (hydro_dimension * wj + xj * wj_dx);
pj->density.wcount += wj;
pj->density.wcount_dh -= xj * wj_dx;
@@ -174,14 +174,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
for (k = 0; k < 3; k++) curlvr[k].v *= ri.v;
rhoi.v = mj.v * wi.v;
- rhoi_dh.v = mj.v * (vec_set1(3.0f) * wi.v + xi.v * wi_dx.v);
+ rhoi_dh.v = mj.v * (vec_set1(hydro_dimension) * wi.v + xi.v * wi_dx.v);
wcounti.v = wi.v;
wcounti_dh.v = xi.v * wi_dx.v;
div_vi.v = mj.v * dvdr.v * wi_dx.v;
for (k = 0; k < 3; k++) rot_vi[k].v = mj.v * curlvr[k].v * wi_dx.v;
rhoj.v = mi.v * wj.v;
- rhoj_dh.v = mi.v * (vec_set1(3.0f) * wj.v + xj.v * wj_dx.v);
+ rhoj_dh.v = mi.v * (vec_set1(hydro_dimension) * wj.v + xj.v * wj_dx.v);
wcountj.v = wj.v;
wcountj_dh.v = xj.v * wj_dx.v;
div_vj.v = mi.v * dvdr.v * wj_dx.v;
@@ -250,7 +250,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
kernel_deval(xi, &wi, &wi_dx);
pi->rho += mj * wi;
- pi->rho_dh -= mj * (3.0f * wi + xi * wi_dx);
+ pi->rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= xi * wi_dx;
@@ -327,7 +327,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
for (k = 0; k < 3; k++) curlvr[k].v *= ri.v;
rhoi.v = mj.v * wi.v;
- rhoi_dh.v = mj.v * (vec_set1(3.0f) * wi.v + xi.v * wi_dx.v);
+ rhoi_dh.v = mj.v * (vec_set1(hydro_dimension) * wi.v + xi.v * wi_dx.v);
wcounti.v = wi.v;
wcounti_dh.v = xi.v * wi_dx.v;
div_vi.v = mj.v * dvdr.v * wi_dx.v;
@@ -358,8 +358,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
float r = sqrtf(r2), ri = 1.0f / r;
float xi, xj;
- float hi_inv, hi2_inv;
- float hj_inv, hj2_inv;
+ float hi_inv, hid_inv;
+ float hj_inv, hjd_inv;
float wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, w, dvdr;
float mi, mj, POrho2i, POrho2j, rhoi, rhoj;
float v_sig, omega_ij, Pi_ij, alpha_ij, tc, v_sig_u;
@@ -376,17 +376,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Get the kernel for hi. */
hi_inv = 1.0f / hi;
- hi2_inv = hi_inv * hi_inv;
+ hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
- wi_dr = hi2_inv * hi2_inv * wi_dx;
+ wi_dr = hid_inv * wi_dx;
/* Get the kernel for hj. */
hj_inv = 1.0f / hj;
- hj2_inv = hj_inv * hj_inv;
+ hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
- wj_dr = hj2_inv * hj2_inv * wj_dx;
+ wj_dr = hjd_inv * wj_dx;
/* Compute dv dot r. */
dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
@@ -457,7 +457,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
vector r, r2, ri;
vector xi, xj;
vector hi, hj, hi_inv, hj_inv;
- vector hi2_inv, hj2_inv;
+ vector hid_inv, hjd_inv;
vector wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, dvdr;
vector w;
vector piPOrho2, pjPOrho2, pirho, pjrho, piu, pju;
@@ -564,19 +564,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
hi.v = vec_load(Hi);
hi_inv.v = vec_rcp(hi.v);
hi_inv.v = hi_inv.v - hi_inv.v * (hi.v * hi_inv.v - vec_set1(1.0f));
- hi2_inv.v = hi_inv.v * hi_inv.v;
+ hid_inv = pow_dimension_plus_one_vec(hi_inv);
xi.v = r.v * hi_inv.v;
kernel_deval_vec(&xi, &wi, &wi_dx);
- wi_dr.v = hi2_inv.v * hi2_inv.v * wi_dx.v;
+ wi_dr.v = hid_inv.v * wi_dx.v;
/* Get the kernel for hj. */
hj.v = vec_load(Hj);
hj_inv.v = vec_rcp(hj.v);
hj_inv.v = hj_inv.v - hj_inv.v * (hj.v * hj_inv.v - vec_set1(1.0f));
- hj2_inv.v = hj_inv.v * hj_inv.v;
+ hjd_inv = pow_dimension_plus_one_vec(hj_inv);
xj.v = r.v * hj_inv.v;
kernel_deval_vec(&xj, &wj, &wj_dx);
- wj_dr.v = hj2_inv.v * hj2_inv.v * wj_dx.v;
+ wj_dr.v = hjd_inv.v * wj_dx.v;
/* Compute dv dot r. */
dvdr.v = ((vi[0].v - vj[0].v) * dx[0].v) + ((vi[1].v - vj[1].v) * dx[1].v) +
@@ -659,8 +659,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
float r = sqrtf(r2), ri = 1.0f / r;
float xi, xj;
- float hi_inv, hi2_inv;
- float hj_inv, hj2_inv;
+ float hi_inv, hid_inv;
+ float hj_inv, hjd_inv;
float wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, w, dvdr;
float /*mi,*/ mj, POrho2i, POrho2j, rhoi, rhoj;
float v_sig, omega_ij, Pi_ij, alpha_ij, tc, v_sig_u;
@@ -677,17 +677,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Get the kernel for hi. */
hi_inv = 1.0f / hi;
- hi2_inv = hi_inv * hi_inv;
+ hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
- wi_dr = hi2_inv * hi2_inv * wi_dx;
+ wi_dr = hid_inv * wi_dx;
/* Get the kernel for hj. */
hj_inv = 1.0f / hj;
- hj2_inv = hj_inv * hj_inv;
+ hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
- wj_dr = hj2_inv * hj2_inv * wj_dx;
+ wj_dr = hjd_inv * wj_dx;
/* Compute dv dot r. */
dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
@@ -752,7 +752,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
vector r, r2, ri;
vector xi, xj;
vector hi, hj, hi_inv, hj_inv;
- vector hi2_inv, hj2_inv;
+ vector hid_inv, hjd_inv;
vector wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, dvdr;
vector w;
vector piPOrho2, pjPOrho2, pirho, pjrho, piu, pju;
@@ -856,19 +856,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
hi.v = vec_load(Hi);
hi_inv.v = vec_rcp(hi.v);
hi_inv.v = hi_inv.v - hi_inv.v * (hi.v * hi_inv.v - vec_set1(1.0f));
- hi2_inv.v = hi_inv.v * hi_inv.v;
+ hid_inv = pow_dimension_plus_one_vec(hi_inv);
xi.v = r.v * hi_inv.v;
kernel_deval_vec(&xi, &wi, &wi_dx);
- wi_dr.v = hi2_inv.v * hi2_inv.v * wi_dx.v;
+ wi_dr.v = hid_inv.v * wi_dx.v;
/* Get the kernel for hj. */
hj.v = vec_load(Hj);
hj_inv.v = vec_rcp(hj.v);
hj_inv.v = hj_inv.v - hj_inv.v * (hj.v * hj_inv.v - vec_set1(1.0f));
- hj2_inv.v = hj_inv.v * hj_inv.v;
+ hjd_inv = pow_dimension_plus_one_vec(hj_inv);
xj.v = r.v * hj_inv.v;
kernel_deval_vec(&xj, &wj, &wj_dx);
- wj_dr.v = hj2_inv.v * hj2_inv.v * wj_dx.v;
+ wj_dr.v = hjd_inv.v * wj_dx.v;
/* Compute dv dot r. */
dvdr.v = ((vi[0].v - vj[0].v) * dx[0].v) + ((vi[1].v - vj[1].v) * dx[1].v) +
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index a1448d8686a30de88e24629558f01661daca223a..60f07c2fc31fd5f38d2929679d0e13beb1cc9131 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -18,7 +18,10 @@
******************************************************************************/
#include "adiabatic_index.h"
+#include "dimension.h"
#include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "kernel_hydro.h"
/**
* @brief Returns the internal energy of a particle
@@ -133,7 +136,14 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
- struct part *restrict p, struct xpart *restrict xp) {}
+ struct part *restrict p, struct xpart *restrict xp) {
+
+ p->ti_begin = 0;
+ p->ti_end = 0;
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
+}
/**
* @brief Prepares a particle for the density calculation.
@@ -169,33 +179,33 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Some smoothing length multiples. */
const float h = p->h;
- const float ih = 1.0f / h;
- const float ih2 = ih * ih;
- const float ih4 = ih2 * ih2;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
/* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root;
- p->rho_dh -= 3.0f * p->mass * kernel_root;
+ p->rho_dh -= hydro_dimension * p->mass * kernel_root;
p->density.wcount += kernel_root;
/* Finish the calculation by inserting the missing h-factors */
- p->rho *= ih * ih2;
- p->rho_dh *= ih4;
+ p->rho *= h_inv_dim;
+ p->rho_dh *= h_inv_dim_plus_one;
p->density.wcount *= kernel_norm;
- p->density.wcount_dh *= ih * kernel_gamma * kernel_norm;
+ p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
const float irho = 1.f / p->rho;
/* Compute the derivative term */
- p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
+ p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
/* Finish calculation of the velocity curl components */
- p->density.rot_v[0] *= ih4 * irho;
- p->density.rot_v[1] *= ih4 * irho;
- p->density.rot_v[2] *= ih4 * irho;
+ p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
+ p->density.rot_v[1] *= h_inv_dim_plus_one * irho;
+ p->density.rot_v[2] *= h_inv_dim_plus_one * irho;
/* Finish calculation of the velocity divergence */
- p->density.div_v *= ih4 * irho;
+ p->density.div_v *= h_inv_dim_plus_one * irho;
}
/**
@@ -308,7 +318,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part *restrict p) {
- p->force.h_dt *= p->h * 0.333333333f;
+ p->force.h_dt *= p->h * hydro_dimension_inv;
p->entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
}
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 323f066db6a047db7ae3401a3525afa941bfc047..d37ac491fe5d8ecdf127c217ca025080daf4bbfd 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -58,7 +58,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the density */
pi->rho += mj * wi;
- pi->rho_dh -= mj * (3.f * wi + ui * wi_dx);
+ pi->rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
/* Compute contribution to the number of neighbours */
pi->density.wcount += wi;
@@ -71,7 +71,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the density */
pj->rho += mi * wj;
- pj->rho_dh -= mi * (3.f * wj + uj * wj_dx);
+ pj->rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
/* Compute contribution to the number of neighbours */
pj->density.wcount += wj;
@@ -191,7 +191,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
/* Compute density of pi. */
rhoi.v = mj.v * wi.v;
- rhoi_dh.v = mj.v * (vec_set1(3.0f) * wi.v + xi.v * wi_dx.v);
+ rhoi_dh.v = mj.v * (vec_set1(hydro_dimension) * wi.v + xi.v * wi_dx.v);
wcounti.v = wi.v;
wcounti_dh.v = xi.v * wi_dx.v;
div_vi.v = mj.v * dvdr.v * wi_dx.v;
@@ -199,7 +199,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
/* Compute density of pj. */
rhoj.v = mi.v * wj.v;
- rhoj_dh.v = mi.v * (vec_set1(3.0f) * wj.v + xj.v * wj_dx.v);
+ rhoj_dh.v = mi.v * (vec_set1(hydro_dimension) * wj.v + xj.v * wj_dx.v);
wcountj.v = wj.v;
wcountj_dh.v = xj.v * wj_dx.v;
div_vj.v = mi.v * dvdr.v * wj_dx.v;
@@ -253,7 +253,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Compute contribution to the density */
pi->rho += mj * wi;
- pi->rho_dh -= mj * (3.f * wi + u * wi_dx);
+ pi->rho_dh -= mj * (hydro_dimension * wi + u * wi_dx);
/* Compute contribution to the number of neighbours */
pi->density.wcount += wi;
@@ -356,7 +356,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
/* Compute density of pi. */
rhoi.v = mj.v * wi.v;
- rhoi_dh.v = mj.v * (vec_set1(3.0f) * wi.v + xi.v * wi_dx.v);
+ rhoi_dh.v = mj.v * (vec_set1(hydro_dimension) * wi.v + xi.v * wi_dx.v);
wcounti.v = wi.v;
wcounti_dh.v = xi.v * wi_dx.v;
div_vi.v = mj.v * dvdr.v * wi_dx.v;
@@ -402,17 +402,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
- const float hi2_inv = hi_inv * hi_inv;
+ const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
const float ui = r * hi_inv;
kernel_deval(ui, &wi, &wi_dx);
- const float wi_dr = hi2_inv * hi2_inv * wi_dx;
+ const float wi_dr = hid_inv * wi_dx;
/* Get the kernel for hj. */
const float hj_inv = 1.0f / hj;
- const float hj2_inv = hj_inv * hj_inv;
+ const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
- const float wj_dr = hj2_inv * hj2_inv * wj_dx;
+ const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
const float P_over_rho2_i = pi->force.P_over_rho2;
@@ -485,7 +485,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
vector r, r2, ri;
vector xi, xj;
vector hi, hj, hi_inv, hj_inv;
- vector hi2_inv, hj2_inv;
+ vector hid_inv, hjd_inv;
vector wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, dvdr;
vector piPOrho, pjPOrho, pirho, pjrho;
vector mi, mj;
@@ -580,19 +580,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
hi.v = vec_load(Hi);
hi_inv.v = vec_rcp(hi.v);
hi_inv.v = hi_inv.v - hi_inv.v * (hi.v * hi_inv.v - vec_set1(1.0f));
- hi2_inv.v = hi_inv.v * hi_inv.v;
+ hid_inv = pow_dimension_plus_one_vec(hi_inv); /* 1/h^(d+1) */
xi.v = r.v * hi_inv.v;
kernel_deval_vec(&xi, &wi, &wi_dx);
- wi_dr.v = hi2_inv.v * hi2_inv.v * wi_dx.v;
+ wi_dr.v = hid_inv.v * wi_dx.v;
/* Get the kernel for hj. */
hj.v = vec_load(Hj);
hj_inv.v = vec_rcp(hj.v);
hj_inv.v = hj_inv.v - hj_inv.v * (hj.v * hj_inv.v - vec_set1(1.0f));
- hj2_inv.v = hj_inv.v * hj_inv.v;
+ hjd_inv = pow_dimension_plus_one_vec(hj_inv); /* 1/h^(d+1) */
xj.v = r.v * hj_inv.v;
kernel_deval_vec(&xj, &wj, &wj_dx);
- wj_dr.v = hj2_inv.v * hj2_inv.v * wj_dx.v;
+ wj_dr.v = hjd_inv.v * wj_dx.v;
/* Compute dv dot r. */
dvdr.v = ((vi[0].v - vj[0].v) * dx[0].v) + ((vi[1].v - vj[1].v) * dx[1].v) +
@@ -677,17 +677,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
- const float hi2_inv = hi_inv * hi_inv;
+ const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
const float ui = r * hi_inv;
kernel_deval(ui, &wi, &wi_dx);
- const float wi_dr = hi2_inv * hi2_inv * wi_dx;
+ const float wi_dr = hid_inv * wi_dx;
/* Get the kernel for hj. */
const float hj_inv = 1.0f / hj;
- const float hj2_inv = hj_inv * hj_inv;
+ const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
- const float wj_dr = hj2_inv * hj2_inv * wj_dx;
+ const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
const float P_over_rho2_i = pi->force.P_over_rho2;
@@ -753,7 +753,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
vector r, r2, ri;
vector xi, xj;
vector hi, hj, hi_inv, hj_inv;
- vector hi2_inv, hj2_inv;
+ vector hid_inv, hjd_inv;
vector wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, dvdr;
vector piPOrho, pjPOrho, pirho, pjrho;
vector mj;
@@ -845,19 +845,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
hi.v = vec_load(Hi);
hi_inv.v = vec_rcp(hi.v);
hi_inv.v = hi_inv.v - hi_inv.v * (hi.v * hi_inv.v - vec_set1(1.0f));
- hi2_inv.v = hi_inv.v * hi_inv.v;
+ hid_inv = pow_dimension_plus_one_vec(hi_inv);
xi.v = r.v * hi_inv.v;
kernel_deval_vec(&xi, &wi, &wi_dx);
- wi_dr.v = hi2_inv.v * hi2_inv.v * wi_dx.v;
+ wi_dr.v = hid_inv.v * wi_dx.v;
/* Get the kernel for hj. */
hj.v = vec_load(Hj);
hj_inv.v = vec_rcp(hj.v);
hj_inv.v = hj_inv.v - hj_inv.v * (hj.v * hj_inv.v - vec_set1(1.0f));
- hj2_inv.v = hj_inv.v * hj_inv.v;
+ hjd_inv = pow_dimension_plus_one_vec(hj_inv);
xj.v = r.v * hj_inv.v;
kernel_deval_vec(&xj, &wj, &wj_dx);
- wj_dr.v = hj2_inv.v * hj2_inv.v * wj_dx.v;
+ wj_dr.v = hjd_inv.v * wj_dx.v;
/* Compute dv dot r. */
dvdr.v = ((vi[0].v - vj[0].v) * dx[0].v) + ((vi[1].v - vj[1].v) * dx[1].v) +
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 1614a483ac3007d37770e04d40833f78d19a8cd2..b89b4f70aea1c29e2df4df3615e34552e4ea3f8d 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -18,6 +18,7 @@
******************************************************************************/
#include "adiabatic_index.h"
+#include "hydro.h"
#include "io_properties.h"
#include "kernel_hydro.h"
@@ -54,8 +55,12 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
float convert_u(struct engine* e, struct part* p) {
- return p->entropy * pow_gamma_minus_one(p->rho) *
- hydro_one_over_gamma_minus_one;
+ return hydro_get_internal_energy(p, 0);
+}
+
+float convert_P(struct engine* e, struct part* p) {
+
+ return hydro_get_pressure(p, 0);
}
/**
@@ -68,7 +73,7 @@ float convert_u(struct engine* e, struct part* p) {
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
- *num_fields = 9;
+ *num_fields = 10;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -90,6 +95,8 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
list[8] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
parts, rho, convert_u);
+ list[9] = io_make_output_field_convert_part(
+ "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, rho, convert_P);
}
/**
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 444920ab2bec940af9c10ea0a6ba471b2ab79762..842a028b9031b3c1cf36e4735089d05e038abdca 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -145,6 +145,11 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part *restrict p, struct xpart *restrict xp) {
+ p->ti_begin = 0;
+ p->ti_end = 0;
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
xp->u_full = p->u;
}
@@ -182,25 +187,25 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Some smoothing length multiples. */
const float h = p->h;
- const float ih = 1.0f / h;
- const float ih2 = ih * ih;
- const float ih4 = ih2 * ih2;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
/* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root;
- p->rho_dh -= 3.0f * p->mass * kernel_root;
+ p->rho_dh -= hydro_dimension * p->mass * kernel_root;
p->density.wcount += kernel_root;
/* Finish the calculation by inserting the missing h-factors */
- p->rho *= ih * ih2;
- p->rho_dh *= ih4;
+ p->rho *= h_inv_dim;
+ p->rho_dh *= h_inv_dim_plus_one;
p->density.wcount *= kernel_norm;
- p->density.wcount_dh *= ih * kernel_gamma * kernel_norm;
+ p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
const float irho = 1.f / p->rho;
/* Compute the derivative term */
- p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
+ p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
}
/**
@@ -281,7 +286,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part *restrict p) {
- p->force.h_dt *= p->h * 0.333333333f;
+ p->force.h_dt *= p->h * hydro_dimension_inv;
}
/**
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 311bdc90744059ea1272c2ca6cf419e5fbc1ba82..9eda45760248ebc0dde18e503abe0b1288659518 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -48,7 +48,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
kernel_deval(xi, &wi, &wi_dx);
pi->rho += mj * wi;
- pi->rho_dh -= mj * (3.f * wi + xi * wi_dx);
+ pi->rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= xi * wi_dx;
@@ -58,7 +58,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
kernel_deval(xj, &wj, &wj_dx);
pj->rho += mi * wj;
- pj->rho_dh -= mi * (3.f * wj + xj * wj_dx);
+ pj->rho_dh -= mi * (hydro_dimension * wj + xj * wj_dx);
pj->density.wcount += wj;
pj->density.wcount_dh -= xj * wj_dx;
}
@@ -93,7 +93,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
kernel_deval(xi, &wi, &wi_dx);
pi->rho += mj * wi;
- pi->rho_dh -= mj * (3.f * wi + xi * wi_dx);
+ pi->rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= xi * wi_dx;
}
@@ -128,19 +128,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
- const float hi2_inv = hi_inv * hi_inv;
+ const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
const float xi = r * hi_inv;
float wi, wi_dx;
kernel_deval(xi, &wi, &wi_dx);
- const float wi_dr = hi2_inv * hi2_inv * wi_dx;
+ const float wi_dr = hid_inv * wi_dx;
/* Get the kernel for hj. */
const float hj_inv = 1.0f / hj;
- const float hj2_inv = hj_inv * hj_inv;
+ const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
const float xj = r * hj_inv;
float wj, wj_dx;
kernel_deval(xj, &wj, &wj_dx);
- const float wj_dr = hj2_inv * hj2_inv * wj_dx;
+ const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
@@ -211,19 +211,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
- const float hi2_inv = hi_inv * hi_inv;
+ const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
const float xi = r * hi_inv;
float wi, wi_dx;
kernel_deval(xi, &wi, &wi_dx);
- const float wi_dr = hi2_inv * hi2_inv * wi_dx;
+ const float wi_dr = hid_inv * wi_dx;
/* Get the kernel for hj. */
const float hj_inv = 1.0f / hj;
- const float hj2_inv = hj_inv * hj_inv;
+ const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
const float xj = r * hj_inv;
float wj, wj_dx;
kernel_deval(xj, &wj, &wj_dx);
- const float wj_dr = hj2_inv * hj2_inv * wj_dx;
+ const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 8dde4a23be3705b4a679c6d17eb8f28c1a8ff6b7..c425ecce6de38b32645b367886e039bb950df68d 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -27,6 +27,7 @@
/* Local headers. */
#include "adiabatic_index.h"
#include "common_io.h"
+#include "dimension.h"
#include "error.h"
#include "hydro.h"
#include "kernel_hydro.h"
@@ -39,8 +40,7 @@ void hydro_props_init(struct hydro_props *p,
/* Kernel properties */
p->eta_neighbours = parser_get_param_float(params, "SPH:resolution_eta");
- const float eta3 = p->eta_neighbours * p->eta_neighbours * p->eta_neighbours;
- p->target_neighbours = eta3 * kernel_norm;
+ p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
p->delta_neighbours = parser_get_param_float(params, "SPH:delta_neighbours");
/* Ghost stuff */
@@ -51,21 +51,26 @@ void hydro_props_init(struct hydro_props *p,
p->CFL_condition = parser_get_param_float(params, "SPH:CFL_condition");
const float max_volume_change = parser_get_opt_param_float(
params, "SPH:max_volume_change", hydro_props_default_volume_change);
- p->log_max_h_change = logf(powf(max_volume_change, 0.33333333333f));
+ p->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
}
void hydro_props_print(const struct hydro_props *p) {
message("Adiabatic index gamma: %f.", hydro_gamma);
- message("Hydrodynamic scheme: %s.", SPH_IMPLEMENTATION);
+
+ message("Hydrodynamic scheme: %s in %dD.", SPH_IMPLEMENTATION,
+ (int)hydro_dimension);
+
message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours (eta=%f).",
kernel_name, p->target_neighbours, p->delta_neighbours,
p->eta_neighbours);
+
message("Hydrodynamic integration: CFL parameter: %.4f.", p->CFL_condition);
+
message(
"Hydrodynamic integration: Max change of volume: %.2f "
"(max|dlog(h)/dt|=%f).",
- powf(expf(p->log_max_h_change), 3.f), p->log_max_h_change);
+ powf(expf(p->log_max_h_change), hydro_dimension), p->log_max_h_change);
if (p->max_smoothing_iterations != hydro_props_default_max_iterations)
message("Maximal iterations in ghost task set to %d (default is %d)",
@@ -76,6 +81,7 @@ void hydro_props_print(const struct hydro_props *p) {
void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
writeAttribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
+ writeAttribute_i(h_grpsph, "Dimension", (int)hydro_dimension);
writeAttribute_s(h_grpsph, "Scheme", SPH_IMPLEMENTATION);
writeAttribute_s(h_grpsph, "Kernel function", kernel_name);
writeAttribute_f(h_grpsph, "Kernel target N_ngb", p->target_neighbours);
diff --git a/src/kernel_hydro.h b/src/kernel_hydro.h
index 7a73bf4f3d51550c5bf9096c0121d0faab819c0b..e8a2b23ac9def1e1938358b2502658660fb91ad5 100644
--- a/src/kernel_hydro.h
+++ b/src/kernel_hydro.h
@@ -20,10 +20,23 @@
#ifndef SWIFT_KERNEL_HYDRO_H
#define SWIFT_KERNEL_HYDRO_H
+/**
+ * @file kernel_hydro.h
+ * @brief Kernel functions for SPH (scalar and vector version).
+ *
+ * All constants and kernel coefficients are taken from table 1 of
+ * Dehnen & Aly, MNRAS, 425, pp. 1062-1082 (2012).
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
#include <math.h>
-/* Includes. */
+/* Local headers. */
#include "const.h"
+#include "dimension.h"
#include "error.h"
#include "inline.h"
#include "vector.h"
@@ -35,8 +48,16 @@
#define kernel_name "Cubic spline (M4)"
#define kernel_degree 3 /* Degree of the polynomial */
#define kernel_ivals 2 /* Number of branches */
+#if defined(HYDRO_DIMENSION_3D)
#define kernel_gamma ((float)(1.825742))
#define kernel_constant ((float)(16. * M_1_PI))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma ((float)(1.778002))
+#define kernel_constant ((float)(80. * M_1_PI / 7.))
+#elif defined(HYDRO_DIMENSION_1D)
+#define kernel_gamma ((float)(1.732051))
+#define kernel_constant ((float)(8. / 3.))
+#endif
static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
__attribute__((aligned(16))) = {3.f, -3.f, 0.f, 0.5f, /* 0 < u < 0.5 */
-1.f, 3.f, -3.f, 1.f, /* 0.5 < u < 1 */
@@ -47,10 +68,18 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
/* Coefficients for the kernel. */
#define kernel_name "Quartic spline (M5)"
-#define kernel_degree 4
-#define kernel_ivals 5
+#define kernel_degree 4 /* Degree of the polynomial */
+#define kernel_ivals 5 /* Number of branches */
+#if defined(HYDRO_DIMENSION_3D)
#define kernel_gamma ((float)(2.018932))
#define kernel_constant ((float)(15625. * M_1_PI / 512.))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma ((float)(1.977173))
+#define kernel_constant ((float)(46875. * M_1_PI / 2398.))
+#elif defined(HYDRO_DIMENSION_1D)
+#define kernel_gamma ((float)(1.936492))
+#define kernel_constant ((float)(3125. / 768.))
+#endif
static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
__attribute__((aligned(16))) = {
6.f, 0.f, -2.4f, 0.f, 0.368f, /* 0 < u < 0.2 */
@@ -65,10 +94,18 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
/* Coefficients for the kernel. */
#define kernel_name "Quintic spline (M6)"
-#define kernel_degree 5
-#define kernel_ivals 3
+#define kernel_degree 5 /* Degree of the polynomial */
+#define kernel_ivals 3 /* Number of branches */
+#if defined(HYDRO_DIMENSION_3D)
#define kernel_gamma ((float)(2.195775))
#define kernel_constant ((float)(2187. * M_1_PI / 40.))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma ((float)(2.158131))
+#define kernel_constant ((float)(15309. * M_1_PI / 478.))
+#elif defined(HYDRO_DIMENSION_1D)
+#define kernel_gamma ((float)(2.121321))
+#define kernel_constant ((float)(243. / 40.))
+#endif
static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
__attribute__((aligned(16))) = {
-10.f, 10.f, 0.f,
@@ -85,10 +122,17 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
/* Coefficients for the kernel. */
#define kernel_name "Wendland C2"
-#define kernel_degree 5
-#define kernel_ivals 1
+#define kernel_degree 5 /* Degree of the polynomial */
+#define kernel_ivals 1 /* Number of branches */
+#if defined(HYDRO_DIMENSION_3D)
#define kernel_gamma ((float)(1.936492))
#define kernel_constant ((float)(21. * M_1_PI / 2.))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma ((float)(1.897367))
+#define kernel_constant ((float)(7. * M_1_PI))
+#elif defined(HYDRO_DIMENSION_1D)
+#error "Wendland C2 kernel not defined in 1D."
+#endif
static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
__attribute__((aligned(16))) = {
4.f, -15.f, 20.f, -10.f, 0.f, 1.f, /* 0 < u < 1 */
@@ -99,10 +143,17 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
/* Coefficients for the kernel. */
#define kernel_name "Wendland C4"
-#define kernel_degree 8
-#define kernel_ivals 1
+#define kernel_degree 8 /* Degree of the polynomial */
+#define kernel_ivals 1 /* Number of branches */
+#if defined(HYDRO_DIMENSION_3D)
#define kernel_gamma ((float)(2.207940))
#define kernel_constant ((float)(495. * M_1_PI / 32.))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma ((float)(2.171239))
+#define kernel_constant ((float)(9. * M_1_PI))
+#elif defined(HYDRO_DIMENSION_1D)
+#error "Wendland C4 kernel not defined in 1D."
+#endif
static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
__attribute__((aligned(16))) = {
11.666667f, -64.f, 140.f, -149.333333f, 70.f,
@@ -115,10 +166,17 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
/* Coefficients for the kernel. */
#define kernel_name "Wendland C6"
-#define kernel_degree 11
-#define kernel_ivals 1
+#define kernel_degree 11 /* Degree of the polynomial */
+#define kernel_ivals 1 /* Number of branches */
+#if defined(HYDRO_DIMENSION_3D)
#define kernel_gamma ((float)(2.449490))
#define kernel_constant ((float)(1365. * M_1_PI / 64.))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma ((float)(2.415230))
+#define kernel_constant ((float)(78. * M_1_PI / 7.))
+#elif defined(HYDRO_DIMENSION_1D)
+#error "Wendland C6 kernel not defined in 1D."
+#endif
static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
__attribute__((aligned(16))) = {
32.f, -231.f, 704.f, -1155.f, 1056.f, -462.f,
@@ -137,30 +195,50 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
/* Ok, now comes the real deal. */
/* First some powers of gamma = H/h */
+#define kernel_gamma_inv ((float)(1. / kernel_gamma))
#define kernel_gamma2 ((float)(kernel_gamma * kernel_gamma))
-#define kernel_gamma3 ((float)(kernel_gamma * kernel_gamma * kernel_gamma))
-#define kernel_gamma4 \
+
+/* define gamma^d, gamma^(d+1), 1/gamma^d and 1/gamma^(d+1) */
+#if defined(HYDRO_DIMENSION_3D)
+#define kernel_gamma_dim ((float)(kernel_gamma * kernel_gamma * kernel_gamma))
+#define kernel_gamma_dim_plus_one \
((float)(kernel_gamma * kernel_gamma * kernel_gamma * kernel_gamma))
-#define kernel_igamma ((float)(1. / kernel_gamma))
-#define kernel_igamma2 ((float)(kernel_igamma * kernel_igamma))
-#define kernel_igamma3 ((float)(kernel_igamma * kernel_igamma * kernel_igamma))
-#define kernel_igamma4 \
- ((float)(kernel_igamma * kernel_igamma * kernel_igamma * kernel_igamma))
+#define kernel_gamma_inv_dim \
+ ((float)(1. / (kernel_gamma * kernel_gamma * kernel_gamma)))
+#define kernel_gamma_inv_dim_plus_one \
+ ((float)(1. / (kernel_gamma * kernel_gamma * kernel_gamma * kernel_gamma)))
+#elif defined(HYDRO_DIMENSION_2D)
+#define kernel_gamma_dim ((float)(kernel_gamma * kernel_gamma))
+#define kernel_gamma_dim_plus_one \
+ ((float)(kernel_gamma * kernel_gamma * kernel_gamma))
+#define kernel_gamma_inv_dim ((float)(1. / (kernel_gamma * kernel_gamma)))
+#define kernel_gamma_inv_dim_plus_one \
+ ((float)(1. / (kernel_gamma * kernel_gamma * kernel_gamma)))
+#elif defined(HYDRO_DIMENSION_1D)
+#define kernel_gamma_dim ((float)(kernel_gamma))
+#define kernel_gamma_dim_plus_one ((float)(kernel_gamma * kernel_gamma))
+#define kernel_gamma_inv_dim ((float)(1. / (kernel_gamma)))
+#define kernel_gamma_inv_dim_plus_one \
+ ((float)(1. / (kernel_gamma * kernel_gamma)))
+#endif
-/* The number of branches */
+/* The number of branches (floating point conversion) */
#define kernel_ivals_f ((float)(kernel_ivals))
/* Kernel self contribution (i.e. W(0,h)) */
-#define kernel_root \
- ((float)(kernel_coeffs[kernel_degree]) * kernel_constant * kernel_igamma3)
+#define kernel_root \
+ ((float)(kernel_coeffs[kernel_degree]) * kernel_constant * \
+ kernel_gamma_inv_dim)
/* Kernel normalisation constant (volume term) */
-#define kernel_norm ((float)(4.0 * M_PI * kernel_gamma3 / 3.0))
+#define kernel_norm ((float)(hydro_dimention_unit_sphere * kernel_gamma_dim))
+
+/* ------------------------------------------------------------------------- */
/**
* @brief Computes the kernel function and its derivative.
*
- * Return 0 if $u > \\gamma = H/h$
+ * Returns garbage if $u > \\gamma = H/h$
*
* @param u The ratio of the distance to the smoothing length $u = x/h$.
* @param W (return) The value of the kernel function $W(x,h)$.
@@ -170,17 +248,17 @@ __attribute__((always_inline)) INLINE static void kernel_deval(
float u, float *restrict W, float *restrict dW_dx) {
/* Go to the range [0,1[ from [0,H[ */
- const float x = u * kernel_igamma;
+ const float x = u * kernel_gamma_inv;
-#if kernel_ivals == 1
- /* Only one branch in this case */
- const float *const coeffs = &kernel_coeffs[0];
-#else
+ //#if kernel_ivals == 1
+ ///* Only one branch in this case */
+ //const float *const coeffs = &kernel_coeffs[0];
+ //#else
/* Pick the correct branch of the kernel */
const int temp = (int)(x * kernel_ivals_f);
const int ind = temp > kernel_ivals ? kernel_ivals : temp;
const float *const coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
-#endif
+ //#endif
/* First two terms of the polynomial ... */
float w = coeffs[0] * x + coeffs[1];
@@ -193,30 +271,32 @@ __attribute__((always_inline)) INLINE static void kernel_deval(
}
/* Return everything */
- *W = w * kernel_constant * kernel_igamma3;
- *dW_dx = dw_dx * kernel_constant * kernel_igamma4;
+ *W = w * kernel_constant * kernel_gamma_inv_dim;
+ *dW_dx = dw_dx * kernel_constant * kernel_gamma_inv_dim_plus_one;
}
/**
* @brief Computes the kernel function.
*
+ * Returns garbage if $u > \\gamma = H/h$
+ *
* @param u The ratio of the distance to the smoothing length $u = x/h$.
* @param W (return) The value of the kernel function $W(x,h)$.
*/
__attribute__((always_inline)) INLINE static void kernel_eval(
float u, float *restrict W) {
/* Go to the range [0,1[ from [0,H[ */
- const float x = u * kernel_igamma;
+ const float x = u * kernel_gamma_inv;
-#if kernel_ivals == 1
- /* Only one branch in this case */
- const float *const coeffs = &kernel_coeffs[0];
-#else
+ //#if kernel_ivals == 1
+ ///* Only one branch in this case */
+ //const float *const coeffs = &kernel_coeffs[0];
+ //#else
/* Pick the correct branch of the kernel */
const int temp = (int)(x * kernel_ivals_f);
const int ind = temp > kernel_ivals ? kernel_ivals : temp;
const float *const coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
-#endif
+ //#endif
/* First two terms of the polynomial ... */
float w = coeffs[0] * x + coeffs[1];
@@ -225,20 +305,24 @@ __attribute__((always_inline)) INLINE static void kernel_eval(
for (int k = 2; k <= kernel_degree; k++) w = x * w + coeffs[k];
/* Return everything */
- *W = w * kernel_constant * kernel_igamma3;
+ *W = w * kernel_constant * kernel_gamma_inv_dim;
}
+/* ------------------------------------------------------------------------- */
+
#ifdef WITH_VECTORIZATION
-static const vector kernel_igamma_vec = FILL_VEC((float)kernel_igamma);
+static const vector kernel_gamma_inv_vec = FILL_VEC((float)kernel_gamma_inv);
static const vector kernel_ivals_vec = FILL_VEC((float)kernel_ivals);
static const vector kernel_constant_vec = FILL_VEC((float)kernel_constant);
-static const vector kernel_igamma3_vec = FILL_VEC((float)kernel_igamma3);
+static const vector kernel_gamma_inv_dim_vec =
+ FILL_VEC((float)kernel_gamma_inv_dim);
-static const vector kernel_igamma4_vec = FILL_VEC((float)kernel_igamma4);
+static const vector kernel_gamma_inv_dim_plus_one_vec =
+ FILL_VEC((float)kernel_gamma_inv_dim_plus_one);
/**
* @brief Computes the kernel function and its derivative (Vectorised version).
@@ -254,7 +338,7 @@ __attribute__((always_inline)) INLINE static void kernel_deval_vec(
/* Go to the range [0,1[ from [0,H[ */
vector x;
- x.v = u->v * kernel_igamma_vec.v;
+ x.v = u->v * kernel_gamma_inv_vec.v;
/* Load x and get the interval id. */
vector ind;
@@ -277,8 +361,9 @@ __attribute__((always_inline)) INLINE static void kernel_deval_vec(
}
/* Return everything */
- w->v = w->v * kernel_constant_vec.v * kernel_igamma3_vec.v;
- dw_dx->v = dw_dx->v * kernel_constant_vec.v * kernel_igamma4_vec.v;
+ w->v = w->v * kernel_constant_vec.v * kernel_gamma_inv_dim_vec.v;
+ dw_dx->v =
+ dw_dx->v * kernel_constant_vec.v * kernel_gamma_inv_dim_plus_one_vec.v;
}
#endif
diff --git a/src/space.c b/src/space.c
index 8b7e5f81b539fb02da44dac67a3575356c08c45e..2a6250e1c19eef77b48e9e54dc813312e64aa32d 100644
--- a/src/space.c
+++ b/src/space.c
@@ -41,8 +41,11 @@
/* Local headers. */
#include "atomic.h"
+#include "const.h"
#include "engine.h"
#include "error.h"
+#include "gravity.h"
+#include "hydro.h"
#include "kernel_hydro.h"
#include "lock.h"
#include "minmax.h"
@@ -1416,6 +1419,59 @@ struct cell *space_getcell(struct space *s) {
return c;
}
+/**
+ * @brief Initialises all the particles by setting them into a valid state
+ *
+ * Calls hydro_first_init_part() on all the particles
+ */
+void space_init_parts(struct space *s) {
+
+ const size_t nr_parts = s->nr_parts;
+ struct part *restrict p = s->parts;
+ struct xpart *restrict xp = s->xparts;
+
+ for (size_t i = 0; i < nr_parts; ++i) {
+
+#ifdef HYDRO_DIMENSION_2D
+ p[i].x[2] = 0.f;
+ p[i].v[2] = 0.f;
+#endif
+
+#ifdef HYDRO_DIMENSION_1D
+ p[i].x[1] = p[i].x[2] = 0.f;
+ p[i].v[1] = p[i].v[2] = 0.f;
+#endif
+
+ hydro_first_init_part(&p[i], &xp[i]);
+ }
+}
+
+/**
+ * @brief Initialises all the g-particles by setting them into a valid state
+ *
+ * Calls gravity_first_init_gpart() on all the particles
+ */
+void space_init_gparts(struct space *s) {
+
+ const size_t nr_gparts = s->nr_gparts;
+ struct gpart *restrict gp = s->gparts;
+
+ for (size_t i = 0; i < nr_gparts; ++i) {
+
+#ifdef HYDRO_DIMENSION_2D
+ gp[i].x[2] = 0.f;
+ gp[i].v_full[2] = 0.f;
+#endif
+
+#ifdef HYDRO_DIMENSION_1D
+ gp[i].x[1] = gp[i].x[2] = 0.f;
+ gp[i].v_full[1] = gp[i].v_full[2] = 0.f;
+#endif
+
+ gravity_first_init_gpart(&gp[i]);
+ }
+}
+
/**
* @brief Split the space into cells given the array of particles.
*
@@ -1549,6 +1605,10 @@ void space_init(struct space *s, const struct swift_params *params,
bzero(s->xparts, Npart * sizeof(struct xpart));
}
+ /* Set the particles in a state where they are ready for a run */
+ space_init_parts(s);
+ space_init_gparts(s);
+
/* Init the space lock. */
if (lock_init(&s->lock) != 0) error("Failed to create space spin-lock.");
diff --git a/src/space.h b/src/space.h
index 180759f2af69dd0d77e6e681229353358d554291..6fe3681c85068979a555ff1d78e32ba7577cf3f0 100644
--- a/src/space.h
+++ b/src/space.h
@@ -162,6 +162,8 @@ void space_split(struct space *s, struct cell *cells, int verbose);
void space_do_split(struct space *s, struct cell *c);
void space_do_parts_sort();
void space_do_gparts_sort();
+void space_init_parts(struct space *s);
+void space_init_gparts(struct space *s);
void space_link_cleanup(struct space *s);
void space_clean(struct space *s);
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 4ef8a070eacb956fe07ba1dad80a8768bba5b5f1..cbb3ec87c735164bb64a67d51e1f891a9030d50b 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -247,9 +247,6 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
part->ti_begin = 0;
part->ti_end = 1;
- xpart->v_full[0] = part->v[0];
- xpart->v_full[1] = part->v[1];
- xpart->v_full[2] = part->v[2];
hydro_first_init_part(part, xpart);
++part;
++xpart;
@@ -448,7 +445,7 @@ int main(int argc, char *argv[]) {
message("Hydro implementation: %s", SPH_IMPLEMENTATION);
message("Smoothing length: h = %f", h * size);
message("Kernel: %s", kernel_name);
- message("Neighbour target: N = %f", h * h * h * kernel_norm);
+ message("Neighbour target: N = %f", pow_dimension(h) * kernel_norm);
message("Density target: rho = %f", rho);
message("div_v target: div = %f", vel == 2 ? 3.f : 0.f);
message("curl_v target: curl = [0., 0., %f]", vel == 3 ? -2.f : 0.f);
@@ -461,6 +458,11 @@ int main(int argc, char *argv[]) {
printf("\n");
+#if !defined(HYDRO_DIMENSION_3D)
+ message("test125cells only useful in 3D. Change parameters in const.h !");
+ return 1;
+#endif
+
/* Build the infrastructure */
struct space space;
space.periodic = 0;
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 3a93d468ef2029b0fc8d07e7edfb3e0abbbe481e..3d1d56d9245f9a9aee9e705754c46bb22ba0ef0a 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -329,7 +329,7 @@ int main(int argc, char *argv[]) {
message("Hydro implementation: %s", SPH_IMPLEMENTATION);
message("Smoothing length: h = %f", h * size);
message("Kernel: %s", kernel_name);
- message("Neighbour target: N = %f", h * h * h * kernel_norm);
+ message("Neighbour target: N = %f", pow_dimension(h) * kernel_norm);
message("Density target: rho = %f", rho);
message("div_v target: div = %f", vel == 2 ? 3.f : 0.f);
message("curl_v target: curl = [0., 0., %f]", vel == 3 ? -2.f : 0.f);