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/CONTRIBUTING.md b/CONTRIBUTING.md
index 883a63c034401ec1fceb477fe33c8342f74a87c1..1248883c0aea4d1ecc3cfeaa219b739ee7712de6 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -1,4 +1,12 @@
The SWIFT source code is using a variation of the 'Google' formatting style.
The script 'format.sh' in the root directory applies the clang-format-3.8
tool with our style choices to all the SWIFT C source file. Please apply
-the formatting script to the files before submitting a merge request.
\ No newline at end of file
+the formatting script to the files before submitting a merge request.
+
+The SWIFT code comes with a series of unit tests that are run automatically
+when a push to the master branch occurs. The suite can be run by doing a `make
+check` in the root directory. Please check that the test suite still
+runs with your changes applied before submitting a merge request and add
+relevant unit tests probing the correctness of new modules. An example of how
+to add a test to the suite can be found by considering the tests/testGreeting
+case.
\ No newline at end of file
diff --git a/configure.ac b/configure.ac
index 483937a9ce4b166410ee42e312b24b13551b5d6a..a798372c2f219a8da71f2b79321c140aba23e790 100644
--- a/configure.ac
+++ b/configure.ac
@@ -469,7 +469,7 @@ if test "$enable_warn" != "no"; then
CFLAGS="$CFLAGS -Wall"
;;
intel)
- CFLAGS="$CFLAGS -w2"
+ CFLAGS="$CFLAGS -w2 -Wunused-variable"
;;
*)
AX_CFLAGS_WARN_ALL
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..050bca39a6b7d6f985d630e057a475945471086a
--- /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..81c016f4775a1e5d693e3c2f78c7d1e5dc532701
--- /dev/null
+++ b/examples/GreshoVortex/run.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e glassPlane_128.hdf5 ]
+then
+ echo "Fetching initial glass file for the 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_2D/getGlass.sh b/examples/SedovBlast_2D/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ae3c977064f5e7a408aa249c5fd9089b3c52ecb1
--- /dev/null
+++ b/examples/SedovBlast_2D/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassPlane_128.hdf5
diff --git a/examples/SedovBlast_2D/makeIC.py b/examples/SedovBlast_2D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..a30e727b4dd42c1f058827959cf12e3b4f152181
--- /dev/null
+++ b/examples/SedovBlast_2D/makeIC.py
@@ -0,0 +1,96 @@
+###############################################################################
+ # 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 random
+from numpy import *
+
+# Generates a swift IC file for the Sedov blast test in a periodic cubic box
+
+# Parameters
+gamma = 5./3. # Gas adiabatic index
+rho0 = 1. # Background density
+P0 = 1.e-6 # Background pressure
+E0= 1. # Energy of the explosion
+N_inject = 21 # Number of particles in which to inject energy
+fileName = "sedov.hdf5"
+
+#---------------------------------------------------
+glass = h5py.File("glassPlane_128.hdf5", "r")
+
+# Read particle positions and h from the glass
+pos = glass["/PartType0/Coordinates"][:,:]
+h = glass["/PartType0/SmoothingLength"][:] * 0.3
+
+numPart = size(h)
+vol = 1.
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = zeros(numPart)
+u = zeros(numPart)
+r = zeros(numPart)
+
+r = sqrt((pos[:,0] - 0.5)**2 + (pos[:,1] - 0.5)**2)
+m[:] = rho0 * vol / numPart
+u[:] = P0 / (rho0 * (gamma - 1))
+
+# Make the central particle detonate
+index = argsort(r)
+u[index[0:N_inject]] = E0 / (N_inject * m[0])
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [1., 1., 1.0]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+
+#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/SedovBlast_2D/plotSolution.py b/examples/SedovBlast_2D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..69e4e1232dd5c14b06e8a705f4add391f1f597f0
--- /dev/null
+++ b/examples/SedovBlast_2D/plotSolution.py
@@ -0,0 +1,280 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+ # 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 2D Sedov blast wave.
+# The script works for a given initial box and dumped energy and computes the solution at a later time t.
+
+# Parameters
+rho_0 = 1. # Background Density
+P_0 = 1.e-6 # Background Pressure
+E_0 = 1. # Energy of the explosion
+gas_gamma = 5./3. # Gas polytropic index
+
+
+# ---------------------------------------------------------------
+# 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])
+
+
+# Read the simulation data
+sim = h5py.File("sedov_%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
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Density"][:]
+
+
+# Now, work our the solution....
+
+from scipy.special import gamma as Gamma
+from numpy import *
+
+def calc_a(g,nu=3):
+ """
+ exponents of the polynomials of the sedov solution
+ g - the polytropic gamma
+ nu - the dimension
+ """
+ a = [0]*8
+
+ a[0] = 2.0 / (nu + 2)
+ a[2] = (1-g) / (2*(g-1) + nu)
+ a[3] = nu / (2*(g-1) + nu)
+ a[5] = 2 / (g-2)
+ a[6] = g / (2*(g-1) + nu)
+
+ a[1] = (((nu+2)*g)/(2.0+nu*(g-1.0)) ) * ( (2.0*nu*(2.0-g))/(g*(nu+2.0)**2) - a[2])
+ a[4] = a[1]*(nu+2) / (2-g)
+ a[7] = (2 + nu*(g-1))*a[1]/(nu*(2-g))
+ return a
+
+def calc_beta(v, g, nu=3):
+ """
+ beta values for the sedov solution (coefficients of the polynomials of the similarity variables)
+ v - the similarity variable
+ g - the polytropic gamma
+ nu- the dimension
+ """
+
+ beta = (nu+2) * (g+1) * array((0.25, (g/(g-1))*0.5,
+ -(2 + nu*(g-1))/2.0 / ((nu+2)*(g+1) -2*(2 + nu*(g-1))),
+ -0.5/(g-1)), dtype=float64)
+
+ beta = outer(beta, v)
+
+ beta += (g+1) * array((0.0, -1.0/(g-1),
+ (nu+2) / ((nu+2)*(g+1) -2.0*(2 + nu*(g-1))),
+ 1.0/(g-1)), dtype=float64).reshape((4,1))
+
+ return beta
+
+
+def sedov(t, E0, rho0, g, n=1000, nu=3):
+ """
+ solve the sedov problem
+ t - the time
+ E0 - the initial energy
+ rho0 - the initial density
+ n - number of points (10000)
+ nu - the dimension
+ g - the polytropic gas gamma
+ """
+ # the similarity variable
+ v_min = 2.0 / ((nu + 2) * g)
+ v_max = 4.0 / ((nu + 2) * (g + 1))
+
+ v = v_min + arange(n) * (v_max - v_min) / (n - 1.0)
+
+ a = calc_a(g, nu)
+ beta = calc_beta(v, g=g, nu=nu)
+ lbeta = log(beta)
+
+ r = exp(-a[0] * lbeta[0] - a[2] * lbeta[1] - a[1] * lbeta[2])
+ rho = ((g + 1.0) / (g - 1.0)) * exp(a[3] * lbeta[1] + a[5] * lbeta[3] + a[4] * lbeta[2])
+ p = exp(nu * a[0] * lbeta[0] + (a[5] + 1) * lbeta[3] + (a[4] - 2 * a[1]) * lbeta[2])
+ u = beta[0] * r * 4.0 / ((g + 1) * (nu + 2))
+ p *= 8.0 / ((g + 1) * (nu + 2) * (nu + 2))
+
+ # we have to take extra care at v=v_min, since this can be a special point.
+ # It is not a singularity, however, the gradients of our variables (wrt v) are.
+ # r -> 0, u -> 0, rho -> 0, p-> constant
+
+ u[0] = 0.0; rho[0] = 0.0; r[0] = 0.0; p[0] = p[1]
+
+ # volume of an n-sphere
+ vol = (pi ** (nu / 2.0) / Gamma(nu / 2.0 + 1)) * power(r, nu)
+
+ # note we choose to evaluate the integral in this way because the
+ # volumes of the first few elements (i.e near v=vmin) are shrinking
+ # very slowly, so we dramatically improve the error convergence by
+ # finding the volumes exactly. This is most important for the
+ # pressure integral, as this is on the order of the volume.
+
+ # (dimensionless) energy of the model solution
+ de = rho * u * u * 0.5 + p / (g - 1)
+ # integrate (trapezium rule)
+ q = inner(de[1:] + de[:-1], diff(vol)) * 0.5
+
+ # the factor to convert to this particular problem
+ fac = (q * (t ** nu) * rho0 / E0) ** (-1.0 / (nu + 2))
+
+ # shock speed
+ shock_speed = fac * (2.0 / (nu + 2))
+ rho_s = ((g + 1) / (g - 1)) * rho0
+ r_s = shock_speed * t * (nu + 2) / 2.0
+ p_s = (2.0 * rho0 * shock_speed * shock_speed) / (g + 1)
+ u_s = (2.0 * shock_speed) / (g + 1)
+
+ r *= fac * t
+ u *= fac
+ p *= fac * fac * rho0
+ rho *= rho0
+ return r, p, rho, u, r_s, p_s, rho_s, u_s, shock_speed
+
+
+# The main properties of the solution
+r_s, P_s, rho_s, v_s, r_shock, _, _, _, _ = sedov(time, E_0, rho_0, gas_gamma, 1000, 2)
+
+# Append points for after the shock
+r_s = np.insert(r_s, np.size(r_s), [r_shock, r_shock*1.5])
+rho_s = np.insert(rho_s, np.size(rho_s), [rho_0, rho_0])
+P_s = np.insert(P_s, np.size(P_s), [P_0, P_0])
+v_s = np.insert(v_s, np.size(v_s), [0, 0])
+
+# Additional arrays
+u_s = P_s / (rho_s * (gas_gamma - 1.)) #internal energy
+s_s = P_s / rho_s**gas_gamma # entropic function
+
+
+
+# Plot the interesting quantities
+figure()
+
+# Velocity profile --------------------------------
+subplot(231)
+plot(r, v_r, '.', color='r', ms=1.)
+plot(r_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Radial~velocity}}~v_r$", labelpad=0)
+xlim(0, 1.3 * r_shock)
+ylim(-0.2, 3.8)
+
+# Density profile --------------------------------
+subplot(232)
+plot(r, rho, '.', color='r', ms=1.)
+plot(r_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=2)
+xlim(0, 1.3 * r_shock)
+ylim(-0.2, 5.2)
+
+# Pressure profile --------------------------------
+subplot(233)
+plot(r, P, '.', color='r', ms=1.)
+plot(r_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(0, 1.3 * r_shock)
+ylim(-1, 12.5)
+
+# Internal energy profile -------------------------
+subplot(234)
+plot(r, u, '.', color='r', ms=1.)
+plot(r_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(0, 1.3 * r_shock)
+ylim(-2, 22)
+
+# Entropy profile ---------------------------------
+subplot(235)
+plot(r, S, '.', color='r', ms=1.)
+plot(r_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+xlim(0, 1.3 * r_shock)
+ylim(-5, 50)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Sedov blast with $\\gamma=%.3f$ in 2D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
+text(-0.49, 0.8, "Background $\\rho_0=%.2f$"%(rho_0), fontsize=10)
+text(-0.49, 0.7, "Energy injected $E_0=%.2f$"%(E_0), 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("Sedov.png", dpi=200)
+
+
+
+
diff --git a/examples/SedovBlast_2D/run.sh b/examples/SedovBlast_2D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..b6f6f73db3ff281492d2abc21eb5ded67e528516
--- /dev/null
+++ b/examples/SedovBlast_2D/run.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e glassPlane_128.hdf5 ]
+then
+ echo "Fetching initial glass file for the Sedov blast example..."
+ ./getGlass.sh
+fi
+if [ ! -e sedov.hdf5 ]
+then
+ echo "Generating initial conditions for the Sedov blast example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 1 sedov.yml
+
+# Plot the solution
+python plotSolution.py 5
diff --git a/examples/SedovBlast_2D/sedov.yml b/examples/SedovBlast_2D/sedov.yml
new file mode 100644
index 0000000000000000000000000000000000000000..6f519835d26ff5aa851ffb8999e650815c522cd3
--- /dev/null
+++ b/examples/SedovBlast_2D/sedov.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: 5e-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-4 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sedov # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1e-2 # 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: ./sedov.hdf5 # The file to read
+
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/examples/main.c b/examples/main.c
index 30e27cc55a6ec508dc352f1acc7f2eed71d71bdd..efd40c3da9802ebfab72a51eb841c035c8490eba 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -455,9 +455,6 @@ int main(int argc, char *argv[]) {
fflush(stdout);
}
- /* Write the state of the system before starting time integration. */
- if (!dry_run) engine_dump_snapshot(&e);
-
/* Init the runner history. */
#ifdef HIST
for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
@@ -493,6 +490,9 @@ int main(int argc, char *argv[]) {
/* Initialise the particles */
engine_init_particles(&e, flag_entropy_ICs);
+ /* Write the state of the system before starting time integration. */
+ engine_dump_snapshot(&e);
+
/* Legend */
if (myrank == 0)
printf("# %6s %14s %14s %10s %10s %16s [%s]\n", "Step", "Time", "Time-step",
diff --git a/examples/plot_tasks.py b/examples/plot_tasks.py
index f2d0aa95d1f35f30476e1989349a07be8d9e5b0a..6e71f476a106937f43bd4bd5973af01f65218afe 100755
--- a/examples/plot_tasks.py
+++ b/examples/plot_tasks.py
@@ -56,8 +56,8 @@ pl.rcParams.update(PLOT_PARAMS)
# Tasks and subtypes. Indexed as in tasks.h.
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair", "init", "ghost",
- "drift", "kick", "kick_fixdt", "send", "recv", "grav_pp", "grav_mm",
- "grav_up", "grav_down", "grav_external", "part_sort", "gpart_sort",
+ "drift", "kick", "kick_fixdt", "send", "recv", "grav_gather_m", "grav_fft",
+ "grav_mm", "grav_up", "grav_external", "part_sort", "gpart_sort",
"split_cell", "rewait", "count"]
TASKCOLOURS = {"none": "black",
@@ -73,10 +73,10 @@ TASKCOLOURS = {"none": "black",
"kick_fixdt": "green",
"send": "yellow",
"recv": "magenta",
- "grav_pp": "mediumorchid",
+ "grav_gather_m": "mediumorchid",
+ "grav_fft": "mediumnightblue",
"grav_mm": "mediumturquoise",
"grav_up": "mediumvioletred",
- "grav_down": "mediumnightblue",
"grav_external": "darkred",
"part_sort": "steelblue",
"gpart_sort": "teal" ,
@@ -84,12 +84,13 @@ TASKCOLOURS = {"none": "black",
"rewait": "olive",
"count": "powerblue"}
-SUBTYPES = ["none", "density", "force", "grav", "count"]
+SUBTYPES = ["none", "density", "force", "grav", "tend", "count"]
SUBCOLOURS = {"none": "black",
"density": "red",
"force": "blue",
"grav": "indigo",
+ "tend": "grey"
"count": "purple"}
# Show docs if help is requested.
diff --git a/examples/plot_tasks_MPI.py b/examples/plot_tasks_MPI.py
index 9a92faf9417c9a302831eb8cb2f4471eb672d59c..7550899da2d4a34a5f73b192cbd7c348426786b7 100755
--- a/examples/plot_tasks_MPI.py
+++ b/examples/plot_tasks_MPI.py
@@ -62,8 +62,8 @@ pl.rcParams.update(PLOT_PARAMS)
# Tasks and subtypes. Indexed as in tasks.h.
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair", "init", "ghost",
- "drift", "kick", "kick_fixdt", "send", "recv", "grav_pp", "grav_mm",
- "grav_up", "grav_down", "grav_external", "part_sort", "gpart_sort",
+ "drift", "kick", "kick_fixdt", "send", "recv", "grav_gather_m", "grav_fft",
+ "grav_mm", "grav_up", "grav_external", "part_sort", "gpart_sort",
"split_cell", "rewait", "count"]
TASKCOLOURS = {"none": "black",
@@ -79,10 +79,10 @@ TASKCOLOURS = {"none": "black",
"kick_fixdt": "green",
"send": "yellow",
"recv": "magenta",
- "grav_pp": "mediumorchid",
+ "grav_gather_m": "mediumorchid",
+ "grav_fft": "mediumnightblue",
"grav_mm": "mediumturquoise",
"grav_up": "mediumvioletred",
- "grav_down": "mediumnightblue",
"grav_external": "darkred",
"part_sort": "steelblue",
"gpart_sort": "teal" ,
@@ -90,12 +90,13 @@ TASKCOLOURS = {"none": "black",
"rewait": "olive",
"count": "powerblue"}
-SUBTYPES = ["none", "density", "force", "grav", "count"]
+SUBTYPES = ["none", "density", "force", "grav", "tend", "count"]
SUBCOLOURS = {"none": "black",
"density": "red",
"force": "blue",
"grav": "indigo",
+ "tend": "grey"
"count": "purple"}
# Show docs if help is requested.
diff --git a/src/Makefile.am b/src/Makefile.am
index d750e2444bc172a836c80ec4d64518908164321b..ace9b6ff553b25d305060ba2e40affd7bdb277e9 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -57,6 +57,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 \
+ 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 ba7a1a64d5c97f1f6d276e5969782351762be4b1..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"
@@ -49,6 +55,10 @@
#define hydro_gamma_minus_one 1.f
#define hydro_one_over_gamma_minus_one 1.f
+#else
+
+#error "An adiabatic index needs to be chosen in const.h !"
+
#endif
/**
@@ -120,8 +130,8 @@ __attribute__((always_inline)) INLINE static float pow_minus_gamma_minus_one(
#if defined(HYDRO_GAMMA_5_3)
- const float inv_cbrt = 1.f / cbrtf(x); /* x^(-1/3) */
- return inv_cbrt * inv_cbrt; /* x^(-2/3) */
+ const float cbrt_inv = 1.f / cbrtf(x); /* x^(-1/3) */
+ return cbrt_inv * cbrt_inv; /* x^(-2/3) */
#elif defined(HYDRO_GAMMA_4_3)
diff --git a/src/const.h b/src/const.h
index 736ba0155e03dc0f595ac06edfd0763b38302e9c..6710fe9e52dfdbede468282d84d4441b923a59fe 100644
--- a/src/const.h
+++ b/src/const.h
@@ -36,11 +36,20 @@
/* 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
//#define HYDRO_GAMMA_2_1
+/* Equation of state choice */
+#define EOS_IDEAL_GAS
+//#define EOS_ISOTHERMAL_GAS
+
/* Kernel function to use */
#define CUBIC_SPLINE_KERNEL
//#define QUARTIC_SPLINE_KERNEL
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 92671258bc725806d4a5825b2c638355f5bb093f..71c21ef96394401eef67f9a783fdbef4238b36c4 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2661,17 +2661,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
new file mode 100644
index 0000000000000000000000000000000000000000..b4a36e8a3ef0bcc281d1f939f89fde08ecf00be9
--- /dev/null
+++ b/src/equation_of_state.h
@@ -0,0 +1,204 @@
+/*******************************************************************************
+ * 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_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"
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "const.h"
+#include "debug.h"
+#include "inline.h"
+
+/* ------------------------------------------------------------------------- */
+#if defined(EOS_IDEAL_GAS)
+
+/**
+ * @brief Returns the internal energy given density and entropy
+ *
+ * Computes \f$u = \frac{S\rho^{\gamma-1} }{\gamma - 1}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_entropy(float density, float entropy) {
+
+ return entropy * pow_gamma_minus_one(density) *
+ hydro_one_over_gamma_minus_one;
+}
+
+/**
+ * @brief Returns the pressure given density and entropy
+ *
+ * Computes \f$P = S\rho^\gamma\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
+ float density, float entropy) {
+
+ return entropy * pow_gamma(density);
+}
+
+/**
+ * @brief Returns the sound speed given density and entropy
+ *
+ * Computes \f$c = \sqrt{\gamma S \rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
+ float density, float entropy) {
+
+ return sqrtf(hydro_gamma * pow_gamma_minus_one(density) * entropy);
+}
+
+/**
+ * @brief Returns the entropy given density and internal energy
+ *
+ * Computes \f$S = \frac{(\gamma - 1)u}{\rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_entropy_from_internal_energy(float density, float u) {
+
+ return hydro_gamma_minus_one * u * pow_minus_gamma_minus_one(density);
+}
+
+/**
+ * @brief Returns the pressure given density and internal energy
+ *
+ * Computes \f$P = (\gamma - 1)u\rho\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_pressure_from_internal_energy(float density, float u) {
+
+ return hydro_gamma_minus_one * u * density;
+}
+
+/**
+ * @brief Returns the sound speed given density and internal energy
+ *
+ * Computes \f$c = \sqrt{\gamma (\gamma - 1) u }\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_soundspeed_from_internal_energy(float density, float u) {
+
+ return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
+}
+
+/* ------------------------------------------------------------------------- */
+#elif defined(EOS_ISOTHERMAL_GAS)
+
+/**
+ * @brief Returns the internal energy given density and entropy
+ *
+ * @param density The density
+ * @param entropy The entropy
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_entropy(float density, float entropy) {
+
+ error("Missing definition !");
+ return 0.f;
+}
+
+/**
+ * @brief Returns the pressure given density and entropy
+ *
+ * @param density The density
+ * @param entropy The entropy
+ */
+__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
+ float density, float entropy) {
+
+ error("Missing definition !");
+ return 0.f;
+}
+
+/**
+ * @brief Returns the entropy given density and internal energy
+ *
+ * @param density The density
+ * @param u The internal energy
+ */
+__attribute__((always_inline)) INLINE static float
+gas_entropy_from_internal_energy(float density, float u) {
+
+ error("Missing definition !");
+ return 0.f;
+}
+
+/**
+ * @brief Returns the pressure given density and internal energy
+ *
+ * @param density The density
+ * @param u The internal energy
+ */
+__attribute__((always_inline)) INLINE static float
+gas_pressure_from_internal_energy(float density, float u) {
+
+ error("Missing definition !");
+ return 0.f;
+}
+
+/**
+ * @brief Returns the sound speed given density and internal energy
+ *
+ * @param density The density
+ * @param u The internal energy
+ */
+__attribute__((always_inline)) INLINE static float
+gas_soundspeed_from_internal_energy(float density, float u) {
+
+ error("Missing definition !");
+ return 0.f;
+}
+
+/* ------------------------------------------------------------------------- */
+#else
+
+#error "An Equation of state needs to be chosen in const.h !"
+
+#endif
+
+#endif /* SWIFT_EQUATION_OF_STATE_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 45b21531e610391fed6693daa07b8907bcdd8a5b..51e09a5d943a7f3782e484289faddece2567f15d 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -19,9 +19,89 @@
#include "adiabatic_index.h"
#include "approx_math.h"
+#include "equation_of_state.h"
#include <float.h>
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ const struct part *restrict p, float dt) {
+
+ return p->u;
+}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+ const struct part *restrict p, float dt) {
+
+ return gas_pressure_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+ const struct part *restrict p, float dt) {
+
+ return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+ const struct part *restrict p, float dt) {
+
+ return p->force.soundspeed;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed internal
+ * energy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param u The new internal energy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
+ struct part *restrict p, float u) {
+
+ p->u = u;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed entropy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param S The new entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_entropy(
+ struct part *restrict p, float S) {
+
+ p->u = gas_internal_energy_from_entropy(p->rho, S);
+}
+
/**
* @brief Computes the hydro time-step of a given particle
*
@@ -57,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.
@@ -93,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;
}
/**
@@ -137,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);
@@ -151,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);
@@ -229,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;
}
/**
@@ -253,51 +341,3 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {}
-
-/**
- * @brief Returns the internal energy of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
- const struct part *restrict p, float dt) {
-
- return p->u;
-}
-
-/**
- * @brief Returns the pressure of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_pressure(
- const struct part *restrict p, float dt) {
-
- return p->force.P_over_rho2 * p->rho * p->rho / p->rho_dh;
-}
-
-/**
- * @brief Returns the entropy of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_entropy(
- const struct part *restrict p, float dt) {
-
- return hydro_gamma_minus_one * p->u * pow_minus_gamma_minus_one(p->rho);
-}
-
-/**
- * @brief Returns the sound speed of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
- const struct part *restrict p, float dt) {
-
- return p->force.soundspeed;
-}
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index 5db3604e1717ba2c14d50c531c0537409501226d..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;
@@ -468,7 +468,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
vector pia[3], pja[3];
vector piu_dt, pju_dt;
vector pih_dt, pjh_dt;
- vector ci, cj, v_sig, vi_sig, vj_sig;
+ vector ci, cj, v_sig;
vector omega_ij, Pi_ij, balsara;
vector pialpha, pjalpha, alpha_ij, v_sig_u, tc;
int j, k;
@@ -503,12 +503,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
pj[2]->force.soundspeed, pj[3]->force.soundspeed,
pj[4]->force.soundspeed, pj[5]->force.soundspeed,
pj[6]->force.soundspeed, pj[7]->force.soundspeed);
- vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
- pi[3]->force.v_sig, pi[4]->force.v_sig, pi[5]->force.v_sig,
- pi[6]->force.v_sig, pi[7]->force.v_sig);
- vj_sig.v = vec_set(pj[0]->force.v_sig, pj[1]->force.v_sig, pj[2]->force.v_sig,
- pj[3]->force.v_sig, pj[4]->force.v_sig, pj[5]->force.v_sig,
- pj[6]->force.v_sig, pj[7]->force.v_sig);
for (k = 0; k < 3; k++) {
vi[k].v = vec_set(pi[0]->v[k], pi[1]->v[k], pi[2]->v[k], pi[3]->v[k],
pi[4]->v[k], pi[5]->v[k], pi[6]->v[k], pi[7]->v[k]);
@@ -544,10 +538,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
pi[2]->force.soundspeed, pi[3]->force.soundspeed);
cj.v = vec_set(pj[0]->force.soundspeed, pj[1]->force.soundspeed,
pj[2]->force.soundspeed, pj[3]->force.soundspeed);
- vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
- pi[3]->force.v_sig);
- vj_sig.v = vec_set(pj[0]->force.v_sig, pj[1]->force.v_sig, pj[2]->force.v_sig,
- pj[3]->force.v_sig);
for (k = 0; k < 3; k++) {
vi[k].v = vec_set(pi[0]->v[k], pi[1]->v[k], pi[2]->v[k], pi[3]->v[k]);
vj[k].v = vec_set(pj[0]->v[k], pj[1]->v[k], pj[2]->v[k], pj[3]->v[k]);
@@ -574,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) +
@@ -669,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;
@@ -687,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] +
@@ -762,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;
@@ -773,7 +763,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
vector pia[3];
vector piu_dt;
vector pih_dt;
- vector ci, cj, v_sig, vi_sig, vj_sig;
+ vector ci, cj, v_sig;
vector omega_ij, Pi_ij, balsara;
vector pialpha, pjalpha, alpha_ij, v_sig_u, tc;
int j, k;
@@ -806,12 +796,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
pj[2]->force.soundspeed, pj[3]->force.soundspeed,
pj[4]->force.soundspeed, pj[5]->force.soundspeed,
pj[6]->force.soundspeed, pj[7]->force.soundspeed);
- vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
- pi[3]->force.v_sig, pi[4]->force.v_sig, pi[5]->force.v_sig,
- pi[6]->force.v_sig, pi[7]->force.v_sig);
- vj_sig.v = vec_set(pj[0]->force.v_sig, pj[1]->force.v_sig, pj[2]->force.v_sig,
- pj[3]->force.v_sig, pj[4]->force.v_sig, pj[5]->force.v_sig,
- pj[6]->force.v_sig, pj[7]->force.v_sig);
for (k = 0; k < 3; k++) {
vi[k].v = vec_set(pi[0]->v[k], pi[1]->v[k], pi[2]->v[k], pi[3]->v[k],
pi[4]->v[k], pi[5]->v[k], pi[6]->v[k], pi[7]->v[k]);
@@ -846,10 +830,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
pi[2]->force.soundspeed, pi[3]->force.soundspeed);
cj.v = vec_set(pj[0]->force.soundspeed, pj[1]->force.soundspeed,
pj[2]->force.soundspeed, pj[3]->force.soundspeed);
- vi_sig.v = vec_set(pi[0]->force.v_sig, pi[1]->force.v_sig, pi[2]->force.v_sig,
- pi[3]->force.v_sig);
- vj_sig.v = vec_set(pj[0]->force.v_sig, pj[1]->force.v_sig, pj[2]->force.v_sig,
- pj[3]->force.v_sig);
for (k = 0; k < 3; k++) {
vi[k].v = vec_set(pi[0]->v[k], pi[1]->v[k], pi[2]->v[k], pi[3]->v[k]);
vj[k].v = vec_set(pj[0]->v[k], pj[1]->v[k], pj[2]->v[k], pj[3]->v[k]);
@@ -876,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 5efce17d78ddc22911ca297a1ef9024d31971ba9..60f07c2fc31fd5f38d2929679d0e13beb1cc9131 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -18,6 +18,93 @@
******************************************************************************/
#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
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ const struct part *restrict p, float dt) {
+
+ const float entropy = p->entropy + p->entropy_dt * dt;
+
+ return gas_internal_energy_from_entropy(p->rho, entropy);
+}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+ const struct part *restrict p, float dt) {
+
+ const float entropy = p->entropy + p->entropy_dt * dt;
+
+ return gas_pressure_from_entropy(p->rho, entropy);
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+ const struct part *restrict p, float dt) {
+
+ return p->entropy + p->entropy_dt * dt;
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+ const struct part *restrict p, float dt) {
+
+ return p->force.soundspeed;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed internal
+ * energy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param u The new internal energy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
+ struct part *restrict p, float u) {
+
+ p->entropy = gas_entropy_from_internal_energy(p->rho, u);
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed entropy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param S The new entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_entropy(
+ struct part *restrict p, float S) {
+
+ p->entropy = S;
+}
/**
* @brief Computes the hydro time-step of a given particle
@@ -49,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.
@@ -85,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;
}
/**
@@ -140,8 +234,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the pressure */
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
- const float pressure =
- (p->entropy + p->entropy_dt * half_dt) * pow_gamma(p->rho);
+ const float pressure = hydro_get_pressure(p, half_dt);
const float irho = 1.f / p->rho;
@@ -200,8 +293,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Drift the pressure */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
- const float pressure =
- (p->entropy + p->entropy_dt * dt_entr) * pow_gamma(p->rho);
+ const float pressure = hydro_get_pressure(p, dt_entr);
const float irho = 1.f / p->rho;
@@ -226,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);
}
@@ -265,56 +357,6 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {
- p->entropy =
- hydro_gamma_minus_one * p->entropy * pow_minus_gamma_minus_one(p->rho);
-}
-
-/**
- * @brief Returns the internal energy of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
- const struct part *restrict p, float dt) {
-
- const float entropy = p->entropy + p->entropy_dt * dt;
-
- return entropy * pow_gamma_minus_one(p->rho) * hydro_one_over_gamma_minus_one;
-}
-
-/**
- * @brief Returns the pressure of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_pressure(
- const struct part *restrict p, float dt) {
-
- return p->force.P_over_rho2 * p->rho * p->rho / p->rho_dh;
-}
-
-/**
- * @brief Returns the entropy of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_entropy(
- const struct part *restrict p, float dt) {
-
- return p->entropy + p->entropy_dt * dt;
-}
-
-/**
- * @brief Returns the sound speed of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
- const struct part *restrict p, float dt) {
-
- return p->force.soundspeed;
+ /* We read u in the entropy field. We now get S from u */
+ p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy);
}
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 f1f284d5d2b56af9d20b24c91eaf4a481f94851f..842a028b9031b3c1cf36e4735089d05e038abdca 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -19,6 +19,94 @@
#include "adiabatic_index.h"
#include "approx_math.h"
+#include "equation_of_state.h"
+
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not internal energy, this function computes the internal
+ * energy from the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ const struct part *restrict p, float dt) {
+
+ return p->u;
+}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+ const struct part *restrict p, float dt) {
+
+ return gas_pressure_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * For implementations where the main thermodynamic variable
+ * is not entropy, this function computes the entropy from
+ * the thermodynamic variable.
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+ const struct part *restrict p, float dt) {
+
+ return gas_entropy_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+ const struct part *restrict p, float dt) {
+
+ return gas_soundspeed_from_internal_energy(p->rho, p->u);
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed internal
+ * energy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param u The new internal energy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
+ struct part *restrict p, float u) {
+
+ p->u = u;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed entropy
+ *
+ * This overrides the current state of the particle but does *not* change its
+ * time-derivatives
+ *
+ * @param p The particle
+ * @param S The new entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_entropy(
+ struct part *restrict p, float S) {
+
+ p->u = gas_internal_energy_from_entropy(p->rho, S);
+}
/**
* @brief Computes the hydro time-step of a given particle
@@ -57,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;
}
@@ -94,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);
}
/**
@@ -193,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;
}
/**
@@ -230,55 +323,3 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {}
-
-/**
- * @brief Returns the internal energy of a particle
- *
- * For implementations where the main thermodynamic variable
- * is not internal energy, this function computes the internal
- * energy from the thermodynamic variable.
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
- const struct part *restrict p, float dt) {
-
- return p->u;
-}
-
-/**
- * @brief Returns the pressure of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_pressure(
- const struct part *restrict p, float dt) {
-
- return p->u * p->rho * hydro_gamma_minus_one;
-}
-
-/**
- * @brief Returns the entropy of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_entropy(
- const struct part *restrict p, float dt) {
-
- return hydro_gamma_minus_one * p->u * pow_minus_gamma_minus_one(p->rho);
-}
-
-/**
- * @brief Returns the sound speed of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
- const struct part *restrict p, float dt) {
-
- return sqrt(p->u * hydro_gamma * hydro_gamma_minus_one);
-}
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..a53465571dbea73c1e2460491500bf6561066e85 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 617ca4bdadd28ebdea7df5d7e18c51757ebfae44..57b7279a508f10e3647105562e5974072a0c57f0 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"
@@ -1420,6 +1423,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.
*
@@ -1553,6 +1609,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 24b40c9161ee81859068d5800e916c786d8f1240..90313be8dbe817d65fbd0e6a8c30c156747594b1 100644
--- a/src/space.h
+++ b/src/space.h
@@ -165,6 +165,8 @@ void space_split(struct space *s, struct cell *cells, int verbose);
void space_split_mapper(void *map_data, int num_elements, void *extra_data);
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/src/task.c b/src/task.c
index d2ce8755b683bb2da386d41c322c3807774ddb47..9e48f8abeb9659ea2674f48620b58da977a6775a 100644
--- a/src/task.c
+++ b/src/task.c
@@ -43,6 +43,7 @@
/* Local headers. */
#include "atomic.h"
#include "error.h"
+#include "inline.h"
#include "lock.h"
/* Task type names. */
@@ -58,8 +59,11 @@ const char *subtaskID_names[task_type_count] = {"none", "density", "force",
/**
* @brief Computes the overlap between the parts array of two given cells.
*/
-size_t task_cell_overlap_part(const struct cell *ci, const struct cell *cj) {
+__attribute__((always_inline)) INLINE static size_t task_cell_overlap_part(
+ const struct cell *ci, const struct cell *cj) {
+
if (ci == NULL || cj == NULL) return 0;
+
if (ci->parts <= cj->parts &&
ci->parts + ci->count >= cj->parts + cj->count) {
return cj->count;
@@ -67,14 +71,18 @@ size_t task_cell_overlap_part(const struct cell *ci, const struct cell *cj) {
cj->parts + cj->count >= ci->parts + ci->count) {
return ci->count;
}
+
return 0;
}
/**
* @brief Computes the overlap between the gparts array of two given cells.
*/
-size_t task_cell_overlap_gpart(const struct cell *ci, const struct cell *cj) {
+__attribute__((always_inline)) INLINE static size_t task_cell_overlap_gpart(
+ const struct cell *ci, const struct cell *cj) {
+
if (ci == NULL || cj == NULL) return 0;
+
if (ci->gparts <= cj->gparts &&
ci->gparts + ci->gcount >= cj->gparts + cj->gcount) {
return cj->gcount;
@@ -82,6 +90,7 @@ size_t task_cell_overlap_gpart(const struct cell *ci, const struct cell *cj) {
cj->gparts + cj->gcount >= ci->gparts + ci->gcount) {
return ci->gcount;
}
+
return 0;
}
@@ -90,7 +99,8 @@ size_t task_cell_overlap_gpart(const struct cell *ci, const struct cell *cj) {
*
* @param t The #task.
*/
-enum task_actions task_acts_on(const struct task *t) {
+__attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
+ const struct task *t) {
switch (t->type) {
diff --git a/src/task.h b/src/task.h
index 3f9985a6baf6b53396e8572c468f005ee52ba4f6..db74c296bb7988d12e474368ed9abffd97ede49e 100644
--- a/src/task.h
+++ b/src/task.h
@@ -107,6 +107,5 @@ int task_lock(struct task *t);
void task_print_mask(unsigned int mask);
void task_print_submask(unsigned int submask);
void task_do_rewait(struct task *t);
-enum task_actions task_acts_on(const struct task *t);
#endif /* SWIFT_TASK_H */
diff --git a/src/timestep.h b/src/timestep.h
index f32a699f6338362c3de247e5a5abc419b9547374..569120cf9cf989b633da35529f7693c8f62a1910 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -106,8 +106,9 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
const float new_dt_external = gravity_compute_timestep_external(
e->external_potential, e->physical_constants, p->gpart);
- const float new_dt_self =
- gravity_compute_timestep_self(e->physical_constants, p->gpart);
+ /* const float new_dt_self = */
+ /* gravity_compute_timestep_self(e->physical_constants, p->gpart); */
+ const float new_dt_self = FLT_MAX; // MATTHIEU
new_dt_grav = fminf(new_dt_external, new_dt_self);
}
diff --git a/src/tools.c b/src/tools.c
index b4a37328ff7df9de76adceebb5ed6801fa417851..b64e17849081994e1969d5f8de0636768dcb729a 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -144,10 +144,7 @@ void pairs_single_density(double *dim, long long int pid,
p = parts[k];
printf("pairs_single: part[%i].id == %lli.\n", k, pid);
- p.rho = 0.0;
- p.density.wcount = 0.0;
- // p.icount = 0;
- p.rho_dh = 0.0;
+ hydro_init_part(&p);
/* Loop over all particle pairs. */
for (k = 0; k < N; k++) {
@@ -487,25 +484,24 @@ void density_dump(int N) {
*/
void engine_single_density(double *dim, long long int pid,
struct part *restrict parts, int N, int periodic) {
- int i, k;
double r2, dx[3];
- float fdx[3], ih;
+ float fdx[3];
struct part p;
/* Find "our" part. */
+ int k;
for (k = 0; k < N && parts[k].id != pid; k++)
;
if (k == N) error("Part not found.");
p = parts[k];
/* Clear accumulators. */
- ih = 1.0f / p.h;
hydro_init_part(&p);
/* Loop over all particle pairs (force). */
- for (k = 0; k < N; k++) {
+ for (int k = 0; k < N; k++) {
if (parts[k].id == p.id) continue;
- for (i = 0; i < 3; i++) {
+ for (int i = 0; i < 3; i++) {
dx[i] = p.x[i] - parts[k].x[i];
if (periodic) {
if (dx[i] < -dim[i] / 2)
@@ -522,12 +518,9 @@ void engine_single_density(double *dim, long long int pid,
}
/* Dump the result. */
- p.rho = ih * ih * ih * (p.rho + p.mass * kernel_root);
- p.rho_dh = p.rho_dh * ih * ih * ih * ih;
- p.density.wcount =
- (p.density.wcount + kernel_root) * (4.0f / 3.0 * M_PI * kernel_gamma3);
- message("part %lli (h=%e) has wcount=%e, rho=%e, rho_dh=%e.", p.id, p.h,
- p.density.wcount, p.rho, p.rho_dh);
+ hydro_end_density(&p, 0);
+ message("part %lli (h=%e) has wcount=%e, rho=%e.", p.id, p.h,
+ p.density.wcount, p.rho);
fflush(stdout);
}
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 8aeee4963ce6e3cc19f4bf6a61999715fd371589..cbb3ec87c735164bb64a67d51e1f891a9030d50b 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -178,6 +178,21 @@ void get_solution(const struct cell *main_cell, struct solution_part *solution,
solution[i].a_hydro[2] = -gradP[2] / solution[i].rho;
solution[i].v_sig = 2.f * solution[i].c;
+
+ solution[i].S_dt = 0.f;
+ solution[i].u_dt = -(solution[i].P / solution[i].rho) * solution[i].div_v;
+ }
+}
+
+void reset_particles(struct cell *c, enum velocity_field vel,
+ enum pressure_field press, float size, float density) {
+
+ for (size_t i = 0; i < c->count; ++i) {
+
+ struct part *p = &c->parts[i];
+
+ set_velocity(p, vel, size);
+ set_energy_state(p, press, size, density);
}
}
@@ -232,9 +247,7 @@ 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;
}
@@ -312,15 +325,16 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
hydro_get_soundspeed(&main_cell->parts[pid], 0.f),
main_cell->parts[pid].a_hydro[0], main_cell->parts[pid].a_hydro[1],
main_cell->parts[pid].a_hydro[2], main_cell->parts[pid].force.h_dt,
- main_cell->parts[pid].force.v_sig,
#if defined(GADGET2_SPH)
- main_cell->parts[pid].entropy_dt, 0.f
+ main_cell->parts[pid].force.v_sig, main_cell->parts[pid].entropy_dt,
+ 0.f
#elif defined(DEFAULT_SPH)
- 0.f, main_cell->parts[pid].force.u_dt
+ main_cell->parts[pid].force.v_sig, 0.f,
+ main_cell->parts[pid].force.u_dt
#elif defined(MINIMAL_SPH)
- 0.f, main_cell->parts[pid].u_dt
+ main_cell->parts[pid].force.v_sig, 0.f, main_cell->parts[pid].u_dt
#else
- 0.f, 0.f
+ 0.f, 0.f, 0.f
#endif
);
}
@@ -340,7 +354,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
solution[pid].div_v, solution[pid].S, solution[pid].u,
solution[pid].P, solution[pid].c, solution[pid].a_hydro[0],
solution[pid].a_hydro[1], solution[pid].a_hydro[2],
- solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt, 0.f);
+ solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt,
+ solution[pid].u_dt);
}
}
@@ -368,7 +383,7 @@ int main(int argc, char *argv[]) {
clocks_set_cpufreq(cpufreq);
/* Choke on FP-exceptions */
- //feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+ feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
/* Get some randomness going */
srand(0);
@@ -427,9 +442,10 @@ int main(int argc, char *argv[]) {
/* Help users... */
message("Adiabatic index: ga = %f", hydro_gamma);
+ 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);
@@ -442,18 +458,28 @@ 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;
space.h_max = h;
+ struct phys_const prog_const;
+ prog_const.const_newton_G = 1.f;
+
struct hydro_props hp;
hp.target_neighbours = h * h * h * kernel_norm;
hp.delta_neighbours = 1.;
hp.max_smoothing_iterations = 1;
+ hp.CFL_condition = 0.1;
struct engine engine;
engine.hydro_properties = &hp;
+ engine.physical_constants = &prog_const;
engine.s = &space;
engine.time = 0.1f;
engine.ti_current = 1;
@@ -508,7 +534,7 @@ int main(int argc, char *argv[]) {
for (int j = 0; j < 125; ++j) runner_do_init(&runner, cells[j], 0);
/* Do the density calculation */
-#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
/* Run all the pairs (only once !)*/
for (int i = 0; i < 5; i++) {
@@ -550,7 +576,7 @@ int main(int argc, char *argv[]) {
for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j]);
/* Do the force calculation */
-#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
/* Do the pairs (for the central 27 cells) */
for (int i = 1; i < 4; i++) {
@@ -585,6 +611,9 @@ int main(int argc, char *argv[]) {
/* Output timing */
message("SWIFT calculation took : %15lli ticks.", time / runs);
+ for (int j = 0; j < 125; ++j)
+ reset_particles(cells[j], vel, press, size, rho);
+
/* NOW BRUTE-FORCE CALCULATION */
const ticks tic = getticks();
@@ -593,7 +622,7 @@ int main(int argc, char *argv[]) {
for (int j = 0; j < 125; ++j) runner_do_init(&runner, cells[j], 0);
/* Do the density calculation */
-#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
/* Run all the pairs (only once !)*/
for (int i = 0; i < 5; i++) {
@@ -634,7 +663,7 @@ int main(int argc, char *argv[]) {
for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j]);
/* Do the force calculation */
-#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
/* Do the pairs (for the central 27 cells) */
for (int i = 1; i < 4; i++) {
@@ -650,6 +679,7 @@ int main(int argc, char *argv[]) {
/* And now the self-interaction for the main cell */
self_all_force(&runner, main_cell);
+
#endif
/* Finally, give a gentle kick */
diff --git a/tests/test125cells.sh.in b/tests/test125cells.sh.in
index 74180108cb89b12e769b6f37b993297d7a163fe1..1d3b0db75d70bf2d5047f71b183812702305df75 100755
--- a/tests/test125cells.sh.in
+++ b/tests/test125cells.sh.in
@@ -7,17 +7,13 @@ do
rm -f brute_force_125_standard.dat swift_dopair_125_standard.dat
- echo ./test125cells -n 6 -r 1 -v $v -p $p -f standard
- echo "running"
./test125cells -n 6 -r 1 -v $v -p $p -f standard
- echo "done"
if [ -e brute_force_125_standard.dat ]
then
python @srcdir@/difffloat.py brute_force_125_standard.dat swift_dopair_125_standard.dat @srcdir@/tolerance_125.dat 6
else
- echo exit 1
- # exit 1
+ exit 1
fi
done
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 09981af2b99d050bb2db3395d19396defe3b8153..3d1d56d9245f9a9aee9e705754c46bb22ba0ef0a 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -147,8 +147,6 @@ void clean_up(struct cell *ci) {
*/
void zero_particle_fields(struct cell *c) {
for (size_t pid = 0; pid < c->count; pid++) {
- c->parts[pid].rho = 0.f;
- c->parts[pid].rho_dh = 0.f;
hydro_init_part(&c->parts[pid]);
}
}
@@ -187,7 +185,12 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
main_cell->parts[pid].x[1], main_cell->parts[pid].x[2],
main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
main_cell->parts[pid].v[2], main_cell->parts[pid].rho,
- main_cell->parts[pid].rho_dh, main_cell->parts[pid].density.wcount,
+#if defined(GIZMO_SPH)
+ 0.f,
+#else
+ main_cell->parts[pid].rho_dh,
+#endif
+ main_cell->parts[pid].density.wcount,
main_cell->parts[pid].density.wcount_dh,
#if defined(GADGET2_SPH)
main_cell->parts[pid].density.div_v,
@@ -223,7 +226,12 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
"%13e %13e %13e\n",
cj->parts[pjd].id, cj->parts[pjd].x[0], cj->parts[pjd].x[1],
cj->parts[pjd].x[2], cj->parts[pjd].v[0], cj->parts[pjd].v[1],
- cj->parts[pjd].v[2], cj->parts[pjd].rho, cj->parts[pjd].rho_dh,
+ cj->parts[pjd].v[2], cj->parts[pjd].rho,
+#if defined(GIZMO_SPH)
+ 0.f,
+#else
+ main_cell->parts[pjd].rho_dh,
+#endif
cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
#if defined(GADGET2_SPH)
cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0],
@@ -318,9 +326,10 @@ int main(int argc, char *argv[]) {
/* Help users... */
message("Adiabatic index: ga = %f", hydro_gamma);
+ 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);
@@ -366,7 +375,7 @@ int main(int argc, char *argv[]) {
const ticks tic = getticks();
-#if defined(DEFAULT_SPH) || defined(GADGET2_SPH)
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
/* Run all the pairs */
for (int j = 0; j < 27; ++j)
@@ -402,7 +411,7 @@ int main(int argc, char *argv[]) {
const ticks tic = getticks();
-#if defined(DEFAULT_SPH) || defined(GADGET2_SPH)
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
/* Run all the brute-force pairs */
for (int j = 0; j < 27; ++j)
diff --git a/tests/testFFT.c b/tests/testFFT.c
index ee916443ba2686bd0d49af15c956101cb301d6e1..c4aeb2885c788bd769bda49bdd15ab121dd8e9d4 100644
--- a/tests/testFFT.c
+++ b/tests/testFFT.c
@@ -23,10 +23,7 @@
#ifndef HAVE_FFTW
-int main() {
-
- return 0;
-}
+int main() { return 0; }
#else
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index 55858b17a51eb37e9e6184b8cf4fc853951bfd57..9d8dca79f3f84ea06e42d61390e68467a5f1b415 100644
--- a/tests/testInteractions.c
+++ b/tests/testInteractions.c
@@ -17,6 +17,12 @@
*
******************************************************************************/
+#include "../config.h"
+
+#ifndef WITH_VECTORIZATION
+int main() { return 0; }
+#else
+
#include <fenv.h>
#include <stdio.h>
#include <stdlib.h>
@@ -138,7 +144,7 @@ void write_header(char *fileName) {
fclose(file);
}
-/*
+/**
* @brief Calls the serial and vectorised version of the non-symmetrical density
* interaction.
*
@@ -324,3 +330,5 @@ int main(int argc, char *argv[]) {
return 0;
}
+
+#endif /* WITH_VECTORIZATION */
diff --git a/tests/testKernelGrav.c b/tests/testKernelGrav.c
index 18140da79e74369b450429a49f3037e17cebb890..2733a4a4d041149499d354ef217ae85b1cd35b7f 100644
--- a/tests/testKernelGrav.c
+++ b/tests/testKernelGrav.c
@@ -77,7 +77,7 @@ int main() {
if (fabsf(gadget_w - swift_w) > 2e-7) {
printf("%2d: r= %f h= %f u= %f Wg(r,h)= %f Ws(r,h)= %f\n", k, r, h, u,
- gadget_w, swift_w);
+ gadget_w, swift_w);
printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w);
return 1;
@@ -102,8 +102,8 @@ int main() {
if (fabsf(gadget_w - swift_w) > 2e-7) {
- printf("%2d: r= %f r_lr= %f u= %f Ws(r)= %f Wg(r)= %f\n", k, r, a_smooth, u,
- swift_w, gadget_w);
+ printf("%2d: r= %f r_lr= %f u= %f Ws(r)= %f Wg(r)= %f\n", k, r, a_smooth,
+ u, swift_w, gadget_w);
printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w);
return 1;
diff --git a/tests/testMaths.c b/tests/testMaths.c
index fbf8da8419cda100c6b25e2f75d099a10f7430b3..3d8f9a8f9db0cf01276eff89aa44157008cbddc6 100644
--- a/tests/testMaths.c
+++ b/tests/testMaths.c
@@ -59,9 +59,9 @@ int main() {
error = 1;
}
- if(error) {
+ if (error) {
printf("%2d: x= %f exp(x)= %e approx_exp(x)=%e abs=%e rel=%e\n", i, x,
- exp_correct, exp_approx, abs, rel);
+ exp_correct, exp_approx, abs, rel);
return 1;
}
}
diff --git a/tests/testPair.c b/tests/testPair.c
index 385e9a9e709b95685d86559e34b3b4c5fb3d3884..8afc2250434de7fdfeaa19a36a213dd9ea79d861 100644
--- a/tests/testPair.c
+++ b/tests/testPair.c
@@ -106,8 +106,6 @@ void clean_up(struct cell *ci) {
*/
void zero_particle_fields(struct cell *c) {
for (size_t pid = 0; pid < c->count; pid++) {
- c->parts[pid].rho = 0.f;
- c->parts[pid].rho_dh = 0.f;
hydro_init_part(&c->parts[pid]);
}
}
@@ -133,9 +131,14 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
"%13e %13e %13e\n",
ci->parts[pid].id, ci->parts[pid].x[0], ci->parts[pid].x[1],
ci->parts[pid].x[2], ci->parts[pid].v[0], ci->parts[pid].v[1],
- ci->parts[pid].v[2], ci->parts[pid].rho, ci->parts[pid].rho_dh,
+ ci->parts[pid].v[2], ci->parts[pid].rho,
+#if defined(GIZMO_SPH)
+ 0.f,
+#else
+ cj->parts[pid].rho_dh,
+#endif
ci->parts[pid].density.wcount, ci->parts[pid].density.wcount_dh,
-#ifdef GADGET2_SPH
+#if defined(GADGET2_SPH) || defined(DEFAULT_SPH)
ci->parts[pid].density.div_v, ci->parts[pid].density.rot_v[0],
ci->parts[pid].density.rot_v[1], ci->parts[pid].density.rot_v[2]
#else
@@ -152,9 +155,14 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
"%13e %13e %13e\n",
cj->parts[pjd].id, cj->parts[pjd].x[0], cj->parts[pjd].x[1],
cj->parts[pjd].x[2], cj->parts[pjd].v[0], cj->parts[pjd].v[1],
- cj->parts[pjd].v[2], cj->parts[pjd].rho, cj->parts[pjd].rho_dh,
+ cj->parts[pjd].v[2], cj->parts[pjd].rho,
+#if defined(GIZMO_SPH)
+ 0.f,
+#else
+ cj->parts[pjd].rho_dh,
+#endif
cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
-#ifdef GADGET2_SPH
+#if defined(GADGET2_SPH) || defined(DEFAULT_SPH)
cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0],
cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
#else
diff --git a/tests/testPair.sh.in b/tests/testPair.sh.in
index 086db14fdd7377d92a35420dffe0ea1dcca561d5..bd7051b060c4acab6cf5a164af1914715856849b 100755
--- a/tests/testPair.sh.in
+++ b/tests/testPair.sh.in
@@ -1,8 +1,11 @@
#!/bin/bash
-rm brute_force_standard.dat swift_dopair_standard.dat
+
+echo ""
+
+rm -f brute_force_standard.dat swift_dopair_standard.dat
./testPair -p 6 -r 1 -d 0 -f standard
-python @srcdir@/difffloat.py brute_force_standard.dat swift_dopair_standard.dat @srcdir@/tolerance_pair.dat
+python @srcdir@/difffloat.py brute_force_standard.dat swift_dopair_standard.dat @srcdir@/tolerance_pair_normal.dat
exit $?
diff --git a/tests/testPairPerturbed.sh.in b/tests/testPairPerturbed.sh.in
index f17d7559e0084727c74a53a9582009a1867d0d0e..9f214e25a098448a906f9da307ea569e327cfdea 100755
--- a/tests/testPairPerturbed.sh.in
+++ b/tests/testPairPerturbed.sh.in
@@ -1,8 +1,11 @@
#!/bin/bash
-rm brute_force_perturbed.dat swift_dopair_perturbed.dat
+
+echo ""
+
+rm -f brute_force_perturbed.dat swift_dopair_perturbed.dat
./testPair -p 6 -r 1 -d 0.1 -f perturbed
-python @srcdir@/difffloat.py brute_force_perturbed.dat swift_dopair_perturbed.dat @srcdir@/tolerance_pair.dat
+python @srcdir@/difffloat.py brute_force_perturbed.dat swift_dopair_perturbed.dat @srcdir@/tolerance_pair_perturbed.dat
exit $?
diff --git a/tests/testParser.sh.in b/tests/testParser.sh.in
index 6c62bc57c041576ea6e807a0a31052ee1dffb4db..9ccfa5078eafca46de406698b4e030bf46ba912c 100755
--- a/tests/testParser.sh.in
+++ b/tests/testParser.sh.in
@@ -1,4 +1,4 @@
#!/bin/bash
-rm parser_output.yml
+rm -f parser_output.yml
./testParser @srcdir@/testParserInput.yaml
diff --git a/tests/tolerance_27_perturbed.dat b/tests/tolerance_27_perturbed.dat
index 913b620dcdf1210230e4fadacaf304b6abf155ac..f0417d845f83d171dacb6b66024cf9a5dc41c6f1 100644
--- a/tests/tolerance_27_perturbed.dat
+++ b/tests/tolerance_27_perturbed.dat
@@ -1,3 +1,3 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.2e-6 2e-6 2e-5 2e-3 2e-6 2e-6 2e-6 2e-6
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-5 1e-4 2e-5 1.5e-4 1.5e-4 1.5e-4
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.2e-6 2e-6 2.1e-5 2e-3 2.1e-6 2e-6 2e-6 2e-6
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-5 1e-4 2e-5 4e-4 4e-4 4e-4
diff --git a/tests/tolerance_pair.dat b/tests/tolerance_pair_normal.dat
similarity index 100%
rename from tests/tolerance_pair.dat
rename to tests/tolerance_pair_normal.dat
diff --git a/tests/tolerance_pair_perturbed.dat b/tests/tolerance_pair_perturbed.dat
new file mode 100644
index 0000000000000000000000000000000000000000..ca58ff45995158e031eca6b60eec498aa6c627ef
--- /dev/null
+++ b/tests/tolerance_pair_perturbed.dat
@@ -0,0 +1,3 @@
+# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-5 1e-5 2e-5 3e-2 1e-5 1e-5 1e-5 1e-5
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-3 4e-4 8e-3 2e-2 1e-4 1.6e-4 1.6e-4 1.6e-4