diff --git a/examples/GreshoVortex/makeIC.py b/examples/GreshoVortex/makeIC.py
index 2f5bebc00ce0f86d3f4f3cccd030cfff5f90d51d..75af3a777682cb56dbd87bab9d125845d06fe153 100644
--- a/examples/GreshoVortex/makeIC.py
+++ b/examples/GreshoVortex/makeIC.py
@@ -1,7 +1,6 @@
###############################################################################
# This file is part of SWIFT.
- # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 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
@@ -19,7 +18,6 @@
##############################################################################
import h5py
-import random
from numpy import *
import sys
@@ -51,7 +49,6 @@ 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)
diff --git a/examples/KelvinHelmoltz/kelvinHelmholtz.yml b/examples/KelvinHelmoltz/kelvinHelmholtz.yml
new file mode 100644
index 0000000000000000000000000000000000000000..38dd16880a209b885f7ad9c30c024988f4d8228f
--- /dev/null
+++ b/examples/KelvinHelmoltz/kelvinHelmholtz.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.5 # 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: kelvinHelmholtz # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.25 # 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.01 # 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: ./kelvinHelmholtz.hdf5 # The file to read
diff --git a/examples/KelvinHelmoltz/makeIC.py b/examples/KelvinHelmoltz/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..5c8632dea52ef301c453cfbf21c35923f12e2d5a
--- /dev/null
+++ b/examples/KelvinHelmoltz/makeIC.py
@@ -0,0 +1,153 @@
+###############################################################################
+ # 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 *
+import sys
+
+# Generates a swift IC file for the Kelvin-Helmholtz vortex in a periodic box
+
+# Parameters
+L2 = 128 # Particles along one edge in the low-density region
+gamma = 5./3. # Gas adiabatic index
+P1 = 2.5 # Central region pressure
+P2 = 2.5 # Outskirts pressure
+v1 = 0.5 # Central region velocity
+v2 = -0.5 # Outskirts vlocity
+rho1 = 2 # Central density
+rho2 = 1 # Outskirts density
+omega0 = 0.1
+sigma = 0.05 / sqrt(2)
+fileOutputName = "kelvinHelmholtz.hdf5"
+#---------------------------------------------------
+
+# Start by generating grids of particles at the two densities
+numPart2 = L2 * L2
+L1 = int(sqrt(numPart2 / rho2 * rho1))
+numPart1 = L1 * L1
+
+#print "N2 =", numPart2, "N1 =", numPart1
+#print "L2 =", L2, "L1 = ", L1
+#print "rho2 =", rho2, "rho1 =", (float(L1*L1)) / (float(L2*L2))
+
+coords1 = zeros((numPart1, 3))
+coords2 = zeros((numPart2, 3))
+h1 = ones(numPart1) * 1.2348 / L1
+h2 = ones(numPart2) * 1.2348 / L2
+m1 = zeros(numPart1)
+m2 = zeros(numPart2)
+u1 = zeros(numPart1)
+u2 = zeros(numPart2)
+vel1 = zeros((numPart1, 3))
+vel2 = zeros((numPart2, 3))
+
+# Particles in the central region
+for i in range(L1):
+ for j in range(L1):
+
+ index = i * L1 + j
+
+ x = i / float(L1) + 1. / (2. * L1)
+ y = j / float(L1) + 1. / (2. * L1)
+
+ coords1[index, 0] = x
+ coords1[index, 1] = y
+ u1[index] = P1 / (rho1 * (gamma-1.))
+ vel1[index, 0] = v1
+
+# Particles in the outskirts
+for i in range(L2):
+ for j in range(L2):
+
+ index = i * L2 + j
+
+ x = i / float(L2) + 1. / (2. * L2)
+ y = j / float(L2) + 1. / (2. * L2)
+
+ coords2[index, 0] = x
+ coords2[index, 1] = y
+ u2[index] = P2 / (rho2 * (gamma-1.))
+ vel2[index, 0] = v2
+
+
+# Now concatenate arrays
+where1 = abs(coords1[:,1]-0.5) < 0.25
+where2 = abs(coords2[:,1]-0.5) > 0.25
+
+coords = append(coords1[where1, :], coords2[where2, :], axis=0)
+
+#print L2*(L2/2), L1*(L1/2)
+#print shape(coords), shape(coords1[where1,:]), shape(coords2[where2,:])
+#print shape(coords), shape(logical_not(coords1[where1,:])), shape(logical_not(coords2[where2,:]))
+
+vel = append(vel1[where1, :], vel2[where2, :], axis=0)
+h = append(h1[where1], h2[where2], axis=0)
+m = append(m1[where1], m2[where2], axis=0)
+u = append(u1[where1], u2[where2], axis=0)
+numPart = size(h)
+ids = linspace(1, numPart, numPart)
+m[:] = (0.5 * rho1 + 0.5 * rho2) / float(numPart)
+
+# Velocity perturbation
+vel[:,1] = omega0 * sin(4*pi*coords[:,0]) * (exp(-(coords[:,1]-0.25)**2 / (2 * sigma**2)) + exp(-(coords[:,1]-0.75)**2 / (2 * sigma**2)))
+
+#File
+fileOutput = h5py.File(fileOutputName, 'w')
+
+# Header
+grp = fileOutput.create_group("/Header")
+grp.attrs["BoxSize"] = [1., 1., 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["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 = fileOutput.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#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.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+#Particle group
+grp = fileOutput.create_group("/PartType0")
+ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = coords
+ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = vel
+ds = grp.create_dataset('Masses', (numPart, 1), 'f')
+ds[()] = m.reshape((numPart,1))
+ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
+ds[()] = h.reshape((numPart,1))
+ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
+ds[()] = u.reshape((numPart,1))
+ds = grp.create_dataset('ParticleIDs', (numPart,1), 'L')
+ds[()] = ids.reshape((numPart,1))
+
+fileOutput.close()
+
+
diff --git a/examples/KelvinHelmoltz/plotSolution.py b/examples/KelvinHelmoltz/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..9191f3ac7ec75c61d5fdab5d347c86222f787fab
--- /dev/null
+++ b/examples/KelvinHelmoltz/plotSolution.py
@@ -0,0 +1,159 @@
+###############################################################################
+ # 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
+P1 = 2.5 # Central region pressure
+P2 = 2.5 # Outskirts pressure
+v1 = 0.5 # Central region velocity
+v2 = -0.5 # Outskirts vlocity
+rho1 = 2 # Central density
+rho2 = 1 # Outskirts density
+
+# ---------------------------------------------------------------
+# 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("kelvinHelmholtz_%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"][:,:]
+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)
+scatter(pos[:,0], pos[:,1], c=vel[:,0], cmap="PuBu", edgecolors='face', s=4, vmin=-1., vmax=1.)
+text(0.97, 0.97, "${\\rm{Velocity~along}}~x$", ha="right", va="top", backgroundcolor="w")
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=0)
+xlim(0, 1)
+ylim(0, 1)
+
+# Radial density profile --------------------------------
+subplot(232)
+scatter(pos[:,0], pos[:,1], c=rho, cmap="PuBu", edgecolors='face', s=4, vmin=0.8, vmax=2.2)
+text(0.97, 0.97, "${\\rm{Density}}$", ha="right", va="top", backgroundcolor="w")
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=0)
+xlim(0, 1)
+ylim(0, 1)
+
+# Radial pressure profile --------------------------------
+subplot(233)
+scatter(pos[:,0], pos[:,1], c=P, cmap="PuBu", edgecolors='face', s=4, vmin=1, vmax=4)
+text(0.97, 0.97, "${\\rm{Pressure}}$", ha="right", va="top", backgroundcolor="w")
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=0)
+xlim(0, 1)
+ylim(0, 1)
+
+# Internal energy profile --------------------------------
+subplot(234)
+scatter(pos[:,0], pos[:,1], c=u, cmap="PuBu", edgecolors='face', s=4, vmin=1.5, vmax=5.)
+text(0.97, 0.97, "${\\rm{Internal~energy}}$", ha="right", va="top", backgroundcolor="w")
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=0)
+xlim(0, 1)
+ylim(0, 1)
+
+# Radial entropy profile --------------------------------
+subplot(235)
+scatter(pos[:,0], pos[:,1], c=S, cmap="PuBu", edgecolors='face', s=4, vmin=0.5, vmax=3.)
+text(0.97, 0.97, "${\\rm{Entropy}}$", ha="right", va="top", backgroundcolor="w")
+xlabel("${\\rm{Position}}~x$", labelpad=0)
+ylabel("${\\rm{Position}}~y$", labelpad=0)
+xlim(0, 1)
+ylim(0, 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, "Kelvin-Helmholtz instability at $t=%.2f$"%(time), fontsize=10)
+text(-0.49, 0.8, "Centre:~~~ $(P, \\rho, v) = (%.3f, %.3f, %.3f)$"%(P1, rho1, v1), fontsize=10)
+text(-0.49, 0.7, "Outskirts: $(P, \\rho, v) = (%.3f, %.3f, %.3f)$"%(P2, rho2, v2), 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("KelvinHelmholtz.png", dpi=200)
diff --git a/examples/KelvinHelmoltz/run.sh b/examples/KelvinHelmoltz/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..4899ca89bc7bbbf72d15d6ecd3961c146a9c9821
--- /dev/null
+++ b/examples/KelvinHelmoltz/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e kelvinHelmholtz.hdf5 ]
+then
+ echo "Generating initial conditions for the Kelvin-Helmholtz example..."
+ python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -t 1 kelvinHelmholtz.yml
+
+# Plot the solution
+python plotSolution.py 6
diff --git a/examples/SedovBlast_2D/makeIC.py b/examples/SedovBlast_2D/makeIC.py
index dc2c574bb76387ef3bebd714423e4612e5838e0d..a30e727b4dd42c1f058827959cf12e3b4f152181 100644
--- a/examples/SedovBlast_2D/makeIC.py
+++ b/examples/SedovBlast_2D/makeIC.py
@@ -31,8 +31,6 @@ E0= 1. # Energy of the explosion
N_inject = 21 # Number of particles in which to inject energy
fileName = "sedov.hdf5"
-#L = 101
-
#---------------------------------------------------
glass = h5py.File("glassPlane_128.hdf5", "r")
@@ -50,10 +48,9 @@ m = zeros(numPart)
u = zeros(numPart)
r = zeros(numPart)
-for i in range(numPart):
- r[i] = sqrt((pos[i,0] - 0.5)**2 + (pos[i,1] - 0.5)**2)
- m[i] = rho0 * vol / numPart
- u[i] = P0 / (rho0 * (gamma - 1))
+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)
diff --git a/examples/SedovBlast_2D/plotSolution.py b/examples/SedovBlast_2D/plotSolution.py
index af7c3b525dc25b2b10399f9cf27553999b28ab92..69e4e1232dd5c14b06e8a705f4add391f1f597f0 100644
--- a/examples/SedovBlast_2D/plotSolution.py
+++ b/examples/SedovBlast_2D/plotSolution.py
@@ -18,7 +18,7 @@
#
##############################################################################
-# Computes the analytical solution of the 3D Sedov blast wave.
+# 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
@@ -85,6 +85,8 @@ P = sim["/PartType0/Pressure"][:]
rho = sim["/PartType0/Density"][:]
+# Now, work our the solution....
+
from scipy.special import gamma as Gamma
from numpy import *
@@ -190,6 +192,8 @@ def sedov(t, E0, rho0, g, n=1000, nu=3):
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
@@ -202,6 +206,8 @@ v_s = np.insert(v_s, np.size(v_s), [0, 0])
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()
diff --git a/examples/SedovBlast_2D/sedov.py b/examples/SedovBlast_2D/sedov.py
deleted file mode 100755
index 2439a7fda91831a7c54350a597db5f026e36fee9..0000000000000000000000000000000000000000
--- a/examples/SedovBlast_2D/sedov.py
+++ /dev/null
@@ -1,212 +0,0 @@
-###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- #
- # 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 numpy as np
-import scipy.integrate as integrate
-import scipy.optimize as optimize
-import os
-
-# Calculate the analytical solution of the Sedov-Taylor shock wave for a given
-# number of dimensions, gamma and time. We assume dimensionless units and a
-# setup with constant density 1, pressure and velocity 0. An energy 1 is
-# inserted in the center at t 0.
-#
-# The solution is a self-similar shock wave, which was described in detail by
-# Sedov (1959). We follow his notations and solution method.
-#
-# The position of the shock at time t is given by
-# r2 = (E/rho1)^(1/(2+nu)) * t^(2/(2+nu))
-# the energy E is related to the inserted energy E0 by E0 = alpha*E, with alpha
-# a constant which has to be calculated by imposing energy conservation.
-#
-# The density for a given radius at a certain time is determined by introducing
-# the dimensionless position coordinate lambda = r/r2. The density profile as
-# a function of lambda is constant in time and given by
-# rho = rho1 * R(V, gamma, nu)
-# and
-# V = V(lambda, gamma, nu)
-#
-# The function V(lambda, gamma, nu) is found by solving a differential equation
-# described in detail in Sedov (1959) chapter 4, section 5. Alpha is calculated
-# from the integrals in section 11 of the same chapter.
-#
-# Numerically, the complete solution requires the use of 3 quadratures and 1
-# root solver, which are implemented using the GNU Scientific Library (GSL).
-# Since some quadratures call functions that themselves contain quadratures,
-# the problem is highly non-linear and complex and takes a while to converge.
-# Therefore, we tabulate the alpha values and profile values the first time
-# a given set of gamma and nu values is requested and reuse these tabulated
-# values.
-#
-# Reference:
-# Sedov (1959): Sedov, L., Similitude et dimensions en mecanique (7th ed.;
-# Moscou: Editions Mir) - french translation of the original
-# book from 1959.
-
-# dimensionless variable z = gamma*P/R as a function of V (and gamma and nu)
-# R is a dimensionless density, while P is a dimensionless pressure
-# z is hence a sort of dimensionless sound speed
-# The expression below corresponds to eq. 11.9 in Sedov (1959), chapter 4
-def _z(V, gamma, nu):
- if V == 2./(nu+2.)/gamma:
- return 0.
- else:
- return (gamma-1.)*V*V*(V-2./(2.+nu))/2./(2./(2.+nu)/gamma-V)
-
-# differential equation that needs to be solved to obtain lambda for a given V
-# corresponds to eq. 5.11 in Sedov (1959), chapter 4 (omega = 0)
-def _dlnlambda_dV(V, gamma, nu):
- nom = _z(V, gamma, nu) - (V-2./(nu+2.))*(V-2./(nu+2.))
- denom = V*(V-1.)*(V-2./(nu+2.))+nu*(2./(nu+2.)/gamma-V)*_z(V, gamma, nu)
- return nom/denom
-
-# dimensionless variable lambda = r/r2 as a function of V (and gamma and nu)
-# found by solving differential equation 5.11 in Sedov (1959), chapter 4
-# (omega = 0)
-def _lambda(V, gamma, nu):
- if V == 2./(nu+2.)/gamma:
- return 0.
- else:
- V0 = 4./(nu+2.)/(gamma+1.)
- integral, err = integrate.quad(_dlnlambda_dV, V, V0, (gamma, nu),
- limit = 8000)
- return np.exp(-integral)
-
-# dimensionless variable R = rho/rho1 as a function of V (and gamma and nu)
-# found by inverting eq. 5.12 in Sedov (1959), chapter 4 (omega = 0)
-# the integration constant C is found by inserting the R, V and z values
-# at the shock wave, where lambda is 1. These correspond to eq. 11.8 in Sedov
-# (1959), chapter 4.
-def _R(V, gamma, nu):
- if V == 2./(nu+2.)/gamma:
- return 0.
- else:
- C = 8.*gamma*(gamma-1.)/(nu+2.)/(nu+2.)/(gamma+1.)/(gamma+1.) \
- *((gamma-1.)/(gamma+1.))**(gamma-2.) \
- *(4./(nu+2.)/(gamma+1.)-2./(nu+2.))
- lambda1 = _lambda(V, gamma, nu)
- lambda5 = lambda1**(nu+2)
- return (_z(V, gamma, nu)*(V-2./(nu+2.))*lambda5/C)**(1./(gamma-2.))
-
-# function of which we need to find the zero point to invert lambda(V)
-def _lambda_min_lambda(V, lambdax, gamma, nu):
- return _lambda(V, gamma, nu) - lambdax
-
-# dimensionless variable V = v*t/r as a function of lambda (and gamma and nu)
-# found by inverting the function lamdba(V) which is found by solving
-# differential equation 5.11 in Sedov (1959), chapter 4 (omega = 0)
-# the invertion is done by searching the zero point of the function
-# lambda_min_lambda defined above
-def _V_inv(lambdax, gamma, nu):
- if lambdax == 0.:
- return 2./(2.+nu)/gamma;
- else:
- return optimize.brentq(_lambda_min_lambda, 2./(nu+2.)/gamma,
- 4./(nu+2.)/(gamma+1.), (lambdax, gamma, nu))
-
-# integrand of the first integral in eq. 11.24 in Sedov (1959), chapter 4
-def _integrandum1(lambdax, gamma, nu):
- V = _V_inv(lambdax, gamma, nu)
- if nu == 2:
- return _R(V, gamma, nu)*V**2*lambdax**3
- else:
- return _R(V, gamma, nu)*V**2*lambdax**4
-
-# integrand of the second integral in eq. 11.24 in Sedov (1959), chapter 4
-def _integrandum2(lambdax, gamma, nu):
- V = _V_inv(lambdax, gamma, nu)
- if V == 2./(nu+2.)/gamma:
- P = 0.
- else:
- P = _z(V, gamma, nu)*_R(V, gamma, nu)/gamma
- if nu == 2:
- return P*lambdax**3
- else:
- return P*lambdax**4
-
-# calculate alpha = E0/E
-# this corresponds to eq. 11.24 in Sedov (1959), chapter 4
-def get_alpha(gamma, nu):
- integral1, err1 = integrate.quad(_integrandum1, 0., 1., (gamma, nu))
- integral2, err2 = integrate.quad(_integrandum2, 0., 1., (gamma, nu))
-
- if nu == 2:
- return np.pi*integral1+2.*np.pi/(gamma-1.)*integral2
- else:
- return 2.*np.pi*integral1+4.*np.pi/(gamma-1.)*integral2
-
-# get the analytical solution for the Sedov-Taylor blastwave given an input
-# energy E, adiabatic index gamma, and number of dimensions nu, at time t and
-# with a maximal outer region radius maxr
-def get_analytical_solution(E, gamma, nu, t, maxr = 1.):
- # we check for the existance of a datafile with precomputed alpha and
- # profile values
- # if it does not exist, we calculate it here and write it out
- # calculation of alpha and the profile takes a while...
- lvec = np.zeros(1000)
- Rvec = np.zeros(1000)
- fname = "sedov_profile_gamma_{gamma}_nu_{nu}.dat".format(gamma = gamma,
- nu = nu)
- if os.path.exists(fname):
- file = open(fname, "r")
- lines = file.readlines()
- alpha = float(lines[0])
- for i in range(len(lines)-1):
- data = lines[i+1].split()
- lvec[i] = float(data[0])
- Rvec[i] = float(data[1])
- else:
- alpha = get_alpha(gamma, nu)
- for i in range(1000):
- lvec[i] = (i+1)*0.001
- V = _V_inv(lvec[i], gamma, nu)
- Rvec[i] = _R(V, gamma, nu)
- file = open(fname, "w")
- file.write("#{alpha}\n".format(alpha = alpha))
- for i in range(1000):
- file.write("{l}\t{R}\n".format(l = lvec[i], R = Rvec[i]))
-
- xvec = np.zeros(1002)
- rhovec = np.zeros(1002)
- if nu == 2:
- r2 = (E/alpha)**0.25*np.sqrt(t)
- else:
- r2 = (E/alpha)**0.2*t**0.4
-
- for i in range(1000):
- xvec[i] = lvec[i]*r2
- rhovec[i] = Rvec[i]
- xvec[1000] = 1.001*r2
- xvec[1001] = maxr
- rhovec[1000] = 1.
- rhovec[1001] = 1.
-
- return xvec, rhovec
-
-def main():
- E = 1.
- gamma = 1.66667
- nu = 2
- t = 0.001
- x, rho = get_analytical_solution(E, gamma, nu, t)
- for i in range(len(x)):
- print x[i], rho[i]
-
-if __name__ == "__main__":
- main()
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 60f07c2fc31fd5f38d2929679d0e13beb1cc9131..d2d4450fa12374a8a8dec624c5e54ba3d47b99aa 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -320,7 +320,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
p->force.h_dt *= p->h * hydro_dimension_inv;
- p->entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
+ p->entropy_dt *=
+ 0.5f * hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
}
/**
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index d37ac491fe5d8ecdf127c217ca025080daf4bbfd..6766e98e6a6ecda12372bee7354b1cd4ca090885 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -469,8 +469,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->force.v_sig = fmaxf(pj->force.v_sig, v_sig);
/* Change in entropy */
- pi->entropy_dt += 0.5f * mj * visc_term * dvdr;
- pj->entropy_dt += 0.5f * mi * visc_term * dvdr;
+ pi->entropy_dt += mj * visc_term * dvdr;
+ pj->entropy_dt += mi * visc_term * dvdr;
}
/**
@@ -631,7 +631,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
pjh_dt.v = mi.v * dvdr.v * ri.v / pirho.v * wj_dr.v;
/* Change in entropy */
- entropy_dt.v = vec_set1(0.5f) * visc_term.v * dvdr.v;
+ entropy_dt.v = visc_term.v * dvdr.v;
/* Store the forces back on the particles. */
for (k = 0; k < VEC_SIZE; k++) {
@@ -644,7 +644,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
pi[k]->force.v_sig = fmaxf(pi[k]->force.v_sig, v_sig.f[k]);
pj[k]->force.v_sig = fmaxf(pj[k]->force.v_sig, v_sig.f[k]);
pi[k]->entropy_dt += entropy_dt.f[k] * mj.f[k];
- pj[k]->entropy_dt -= entropy_dt.f[k] * mi.f[k];
+ pj[k]->entropy_dt += entropy_dt.f[k] * mi.f[k];
}
#else
@@ -738,7 +738,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
/* Change in entropy */
- pi->entropy_dt += 0.5f * mj * visc_term * dvdr;
+ pi->entropy_dt += mj * visc_term * dvdr;
}
/**
@@ -894,7 +894,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
pih_dt.v = mj.v * dvdr.v * ri.v / pjrho.v * wi_dr.v;
/* Change in entropy */
- entropy_dt.v = vec_set1(0.5f) * mj.v * visc_term.v * dvdr.v;
+ entropy_dt.v = mj.v * visc_term.v * dvdr.v;
/* Store the forces back on the particles. */
for (k = 0; k < VEC_SIZE; k++) {
diff --git a/tests/difffloat.py b/tests/difffloat.py
index 57707c5920997e3ef688606a0839f59a69d2e4f2..e0f0864372264899c6de1bf2f83ab678b7dd9ead 100644
--- a/tests/difffloat.py
+++ b/tests/difffloat.py
@@ -89,13 +89,12 @@ for i in range(n_lines_to_check):
abs_diff = abs(data1[i,j] - data2[i,j])
sum = abs(data1[i,j] + data2[i,j])
- if abs(data1[i,j]) + abs(data2[i,j]) < 2.5e-7: continue
if sum > 0:
rel_diff = abs(data1[i,j] - data2[i,j]) / sum
else:
rel_diff = 0.
- if( abs_diff > absTol[j]):
+ if( abs_diff > 1.1*absTol[j]):
print "Absolute difference larger than tolerance (%e) for particle %d, column %d:"%(absTol[j], i,j)
print "%10s: a = %e"%("File 1", data1[i,j])
print "%10s: b = %e"%("File 2", data2[i,j])
@@ -103,7 +102,9 @@ for i in range(n_lines_to_check):
print ""
error = True
- if( rel_diff > relTol[j]):
+ if abs(data1[i,j]) < 1e-6 and + abs(data2[i,j]) < 1e-6 : continue
+
+ if( rel_diff > 1.1*relTol[j]):
print "Relative difference larger than tolerance (%e) for particle %d, column %d:"%(relTol[j], i,j)
print "%10s: a = %e"%("File 1", data1[i,j])
print "%10s: b = %e"%("File 2", data2[i,j])
diff --git a/tests/test125cells.c b/tests/test125cells.c
index e75754999ccd11e61a930003592169eb67029a7b..4619d1fce34e67d4a1f62af59792390310dcfccb 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -511,8 +511,8 @@ int main(int argc, char *argv[]) {
prog_const.const_newton_G = 1.f;
struct hydro_props hp;
- hp.target_neighbours = h * h * h * kernel_norm;
- hp.delta_neighbours = 1.;
+ hp.target_neighbours = pow_dimension(h) * kernel_norm;
+ hp.delta_neighbours = 2.;
hp.max_smoothing_iterations = 1;
hp.CFL_condition = 0.1;