Merge remote-tracking branch 'origin/master' into faster_rebuilds

Conflicts:
	src/space.c

INSTALL.swift

@@ -34,13 +34,17 @@ directory. See README for run parameters.
 
 SWIFT has been successfully built and tested with the following compilers:
 
-  - GCC 4.8.x  
+  - GCC 4.8.x
   - Intel ICC 15.0.x
-  - clang 3.4.x 
+  - clang 3.4.x
 
 More recent versions and slightly older ones should also be able to
 build the software.
 
+It has also been built with Intel and GNU C++ compilers, but that currently
+requires the --disable-vec and, for Intel, --disable-compiler-warnings
+configure options.
+
 By default an attempt to choose suitable set of optimizing compiler flags
 will be made, targeted for the host machine of the build. If this doesn't
 work or the binaries will for another architecture then you can stop the
@@ -61,7 +65,7 @@ You could also add some additional flags:
 
     ./configure --enable-debug --disable-optimization CFLAGS="-O2"
 
-for instance. GCC address sanitizer flags can be included using the 
+for instance. GCC address sanitizer flags can be included using the
 
     ./configure --enable-sanitizer
 
@@ -113,7 +117,7 @@ before you can build it.
 	none-standard name.
 
 
- - libtool: The build system relies on libtool.
+ - libtool: The build system relies on libtool as well as the other autotools.
 
 
                            Optional Dependencies

configure.ac

@@ -54,6 +54,9 @@ AC_USE_SYSTEM_EXTENSIONS
 AX_COMPILER_VENDOR
 AX_COMPILER_VERSION
 
+#  Restrict support.
+AC_C_RESTRICT
+
 # Interprocedural optimization support. Needs special handling for linking and
 # archiving as well as compilation with Intels, needs to be done before
 # libtool is configured (to use correct LD).
@@ -454,6 +457,14 @@ AC_SUBST([METIS_LIBS])
 AC_SUBST([METIS_INCS])
 AM_CONDITIONAL([HAVEMETIS],[test -n "$METIS_LIBS"])
 
+# METIS fixed width integer printing can require this, so define. Only needed
+# for some non C99 compilers, i.e. C++ pre C++11.
+AH_VERBATIM([__STDC_FORMAT_MACROS],
+            [/* Needed to get PRIxxx macros from stdint.h when not using C99 */
+#ifndef __STDC_FORMAT_MACROS
+#define __STDC_FORMAT_MACROS 1
+#endif])
+
 # Check for grackle. 
 have_grackle="no"
 AC_ARG_WITH([grackle],

examples/CoolingHalo/cooling_halo.yml

 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.9885e39     # 10^6 solar masses
+  UnitMass_in_cgs:     1.98848e39    # 10^6 solar masses
   UnitLength_in_cgs:   3.0856776e21  # Kiloparsecs
   UnitVelocity_in_cgs: 1e5           # Kilometres per second
   UnitCurrent_in_cgs:  1   # Amperes

examples/CoolingHaloWithSpin/cooling_halo.yml

 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.9885e39     # 10^6 solar masses
+  UnitMass_in_cgs:     1.98848e39    # 10^6 solar masses
   UnitLength_in_cgs:   3.0856776e21  # Kiloparsecs
   UnitVelocity_in_cgs: 1e5           # Kilometres per second
   UnitCurrent_in_cgs:  1   # Amperes

examples/EAGLE_100/README

@@ -13,4 +13,4 @@ The particle load of the main EAGLE simulation can be reproduced by
 running these ICs on 4096 cores.
 
 MD5 checksum of the ICs:
-2301ea73e14207b541bbb04163c5269e  EAGLE_ICs_100.hdf5
+bb5c3531e8573739ff4ee35049750ac9  EAGLE_ICs_100.hdf5

examples/EAGLE_100/eagle_100.yml

@@ -10,7 +10,7 @@ InternalUnitSystem:
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
   time_end:   1e-2  # The end time of the simulation (in internal units).
-  dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
+  dt_min:     1e-12 # 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).
 
 Scheduler:
@@ -40,4 +40,3 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  ./EAGLE_ICs_100.hdf5     # The file to read
-  cleanup_h:  1

examples/EvrardCollapse_3D/evrard.yml

+# 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.8   # The end time of the simulation (in internal units).
+  dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-3  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            evrard # Common part of the name of output files
+  time_first:          0.       # Time of the first output (in internal units)
+  delta_time:          0.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).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                   0.025    # Constant dimensionless multiplier for time integration.
+  epsilon:               0.001      # Softening length (in internal units).
+  theta:                 0.9
+  a_smooth:              1.25     # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
+  r_cut:                 4.5      # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./evrard.hdf5       # The file to read
+
+PhysicalConstants:
+  G: 1.

examples/EvrardCollapse_3D/getReference.sh

+#! /bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/evrardCollapse3D_exact.txt

examples/EvrardCollapse_3D/makeIC.py

+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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 Evrard collapse
+
+# Parameters
+gamma = 5. / 3.      # Gas adiabatic index
+M = 1.  # total mass of the sphere
+R = 1.               # radius of the sphere
+u0 = 0.05 / M        # initial thermal energy
+fileName = "evrard.hdf5" 
+numPart = 100000
+
+r = R * sqrt(random.random(numPart))
+phi = 2. * pi * random.random(numPart)
+cos_theta = 2. * random.random(numPart) - 1.
+
+sin_theta = sqrt(1. - cos_theta**2)
+cos_phi = cos(phi)
+sin_phi = sin(phi)
+
+pos = zeros((numPart, 3))
+pos[:,0] = r * sin_theta * cos_phi
+pos[:,1] = r * sin_theta * sin_phi
+pos[:,2] = r * cos_theta
+
+# shift particles to put the sphere in the centre of the box
+pos += array([50. * R, 50. * R, 50. * R])
+
+h = ones(numPart) * 2. * R / numPart**(1. / 3.)
+
+# Generate extra arrays
+v = zeros((numPart, 3))
+ids = linspace(1, numPart, numPart)
+m = ones(numPart) * M / numPart
+u = ones(numPart) * u0
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [100. * R, 100. * R, 100. * R]
+grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 3
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 0
+
+#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()

examples/EvrardCollapse_3D/plotSolution.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#               2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+# 
+# 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/>.
+# 
+################################################################################
+
+# Compares the swift result for the 2D spherical Sod shock with a high
+# resolution 2D reference result
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+from scipy import stats
+import h5py
+
+# 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
+
+# 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("evrard_%04d.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"]
+
+coords = sim["/PartType0/Coordinates"]
+x = sqrt((coords[:,0] - 0.5 * boxSize)**2 + (coords[:,1] - 0.5 * boxSize)**2 + \
+         (coords[:,2] - 0.5 * boxSize)**2)
+vels = sim["/PartType0/Velocities"]
+v = sqrt(vels[:,0]**2 + vels[:,1]**2 + vels[:,2]**2)
+u = sim["/PartType0/InternalEnergy"][:]
+S = sim["/PartType0/Entropy"][:]
+P = sim["/PartType0/Pressure"][:]
+rho = sim["/PartType0/Density"][:]
+
+# Bin the data
+x_bin_edge = logspace(-3., log10(2.), 100)
+x_bin = 0.5*(x_bin_edge[1:] + x_bin_edge[:-1])
+rho_bin,_,_ = stats.binned_statistic(x, rho, statistic='mean', bins=x_bin_edge)
+v_bin,_,_ = stats.binned_statistic(x, v, statistic='mean', bins=x_bin_edge)
+P_bin,_,_ = stats.binned_statistic(x, P, statistic='mean', bins=x_bin_edge)
+S_bin,_,_ = stats.binned_statistic(x, S, statistic='mean', bins=x_bin_edge)
+u_bin,_,_ = stats.binned_statistic(x, u, statistic='mean', bins=x_bin_edge)
+rho2_bin,_,_ = stats.binned_statistic(x, rho**2, statistic='mean', bins=x_bin_edge)
+v2_bin,_,_ = stats.binned_statistic(x, v**2, statistic='mean', bins=x_bin_edge)
+P2_bin,_,_ = stats.binned_statistic(x, P**2, statistic='mean', bins=x_bin_edge)
+S2_bin,_,_ = stats.binned_statistic(x, S**2, statistic='mean', bins=x_bin_edge)
+u2_bin,_,_ = stats.binned_statistic(x, u**2, statistic='mean', bins=x_bin_edge)
+rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
+v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
+P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
+S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
+u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
+
+ref = loadtxt("evrardCollapse3D_exact.txt")
+
+# Plot the interesting quantities
+figure()
+
+# Velocity profile --------------------------------
+subplot(231)
+semilogx(x, -v, '.', color='r', ms=0.2)
+semilogx(ref[:,0], ref[:,2], "k--", alpha=0.8, lw=1.2)
+errorbar(x_bin, -v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Velocity}}~v_r$", labelpad=0)
+xlim(1.e-3, 2.)
+ylim(-1.7, 0.1)
+
+# Density profile --------------------------------
+subplot(232)
+loglog(x, rho, '.', color='r', ms=0.2)
+loglog(ref[:,0], ref[:,1], "k--", alpha=0.8, lw=1.2)
+errorbar(x_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
+xlim(1.e-3, 2.)
+ylim(1.e-2, 1.e4)
+
+# Pressure profile --------------------------------
+subplot(233)
+loglog(x, P, '.', color='r', ms=0.2)
+loglog(ref[:,0], ref[:,3], "k--", alpha=0.8, lw=1.2)
+errorbar(x_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Pressure}}~P$", labelpad=0)
+xlim(1.e-3, 2.)
+ylim(1.e-4, 1.e3)
+
+# Internal energy profile -------------------------
+subplot(234)
+loglog(x, u, '.', color='r', ms=0.2)
+loglog(ref[:,0], ref[:,3] / ref[:,1] / (gas_gamma - 1.), "k--", alpha=0.8, lw=1.2)
+errorbar(x_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
+xlim(1.e-3, 2.)
+ylim(1.e-2, 2.)
+
+# Entropy profile ---------------------------------
+subplot(235)
+semilogx(x, S, '.', color='r', ms=0.2)
+semilogx(ref[:,0], ref[:,3] / ref[:,1]**gas_gamma, "k--", alpha=0.8, lw=1.2)
+errorbar(x_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+ylabel("${\\rm{Entropy}}~S$", labelpad=0)
+xlim(1.e-3, 2.)
+ylim(0., 0.25)
+
+# Information -------------------------------------
+subplot(236, frameon=False)
+
+text(-0.49, 0.9, "Evrard collapse with $\\gamma=%.3f$ in 3D\nat $t=%.2f$"%(gas_gamma,time), 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([])
+
+tight_layout()
+savefig("EvrardCollapse.png", dpi=200)

examples/EvrardCollapse_3D/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e evrard.hdf5 ]
+then
+    echo "Generating initial conditions for the Evrard collapse example..."
+    python makeIC.py
+fi
+
+# Run SWIFT
+../swift -s -G -t 4 evrard.yml 2>&1 | tee output.log
+
+# Get the high resolution 1D reference result if not present.
+if [ ! -e evrardCollapse3D_exact.txt ]
+then
+    echo "Fetching the reference result for the Evrard collapse example..."
+    ./getReference.sh
+fi
+
+# Plot the solution
+python plotSolution.py 8

examples/ExternalPointMass/externalPointMass.yml

 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.9885e33     # Grams
-  UnitLength_in_cgs:   3.0856776e21  # Centimeters
-  UnitVelocity_in_cgs: 1e5           # Centimeters per second
+  UnitMass_in_cgs:     1.98848e33    # M_sun
+  UnitLength_in_cgs:   3.0856776e21  # kpc
+  UnitVelocity_in_cgs: 1e5           # km/s
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
 

examples/ExternalPointMass/makeIC.py

@@ -28,7 +28,7 @@ import random
 
 # physical constants in cgs
 NEWTON_GRAVITY_CGS  = 6.67408e-8
-SOLAR_MASS_IN_CGS   = 1.9885e33
+SOLAR_MASS_IN_CGS   = 1.98848e33
 PARSEC_IN_CGS       = 3.0856776e18
 
 # choice of units
@@ -86,7 +86,7 @@ grp.attrs["PeriodicBoundariesOn"] = periodic
 #Units
 grp = file.create_group("/Units")
 grp.attrs["Unit length in cgs (U_L)"] = 3.0856776e21
-grp.attrs["Unit mass in cgs (U_M)"] = 1.9885e33
+grp.attrs["Unit mass in cgs (U_M)"] = 1.98848e33
 grp.attrs["Unit time in cgs (U_t)"] = 3.0856776e16
 grp.attrs["Unit current in cgs (U_I)"] = 1.
 grp.attrs["Unit temperature in cgs (U_T)"] = 1.

examples/GreshoVortex_3D/getGlass.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5

examples/GreshoVortex_3D/gresho.yml

+# 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
+
+Scheduler:
+  max_top_level_cells: 15
+
+# 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).
+  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

examples/GreshoVortex_3D/makeIC.py

+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#               2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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 Gresho-Chan vortex in a periodic box
+
+# Parameters
+gamma = 5./3.     # Gas adiabatic index
+rho0 = 1           # Gas density
+P0 = 0.           # Constant additional pressure (should have no impact on the dynamics)
+fileOutputName = "greshoVortex.hdf5" 
+fileGlass = "glassCube_64.hdf5"
+#---------------------------------------------------
+
+# 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**3 / numPart
+u = zeros(numPart)
+v = zeros((numPart, 3))
+
+for i in range(numPart):
+    
+    x = coords[i,0]
+    y = coords[i,1]
+
+    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
+fileOutput = h5py.File(fileOutputName, 'w')
+
+# Header
+grp = fileOutput.create_group("/Header")
+grp.attrs["BoxSize"] = [boxSize, 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["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]
+grp.attrs["Dimension"] = 3
+
+#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[()] = v
+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()

examples/GreshoVortex_3D/plotSolution.py

+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#               2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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 *
+from scipy import stats
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (9.90,6.45),
+'figure.subplot.left'    : 0.045,
+'figure.subplot.right'   : 0.99,
+'figure.subplot.bottom'  : 0.05,
+'figure.subplot.top'     : 0.99,
+'figure.subplot.wspace'  : 0.15,
+'figure.subplot.hspace'  : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+snap = int(sys.argv[1])
+
+# 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_%04d.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"][:]
+
+# Bin te data
+r_bin_edge = np.arange(0., 1., 0.02)
+r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
+rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
+v_bin,_,_ = stats.binned_statistic(r, v_phi, statistic='mean', bins=r_bin_edge)
+P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
+S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
+u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
+rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
+v2_bin,_,_ = stats.binned_statistic(r, v_phi**2, statistic='mean', bins=r_bin_edge)
+P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
+S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
+u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_bin_edge)
+rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
+v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
+P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
+S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
+u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
+
+
+# Plot the interesting quantities
+figure()
+
+
+# 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)
+errorbar(r_bin, v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
+plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+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)
+errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
+plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+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)
+errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
+plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+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)
+errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
+plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+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)
+errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
+plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
+plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
+xlabel("${\\rm{Radius}}~r$", labelpad=0)
+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)

examples/GreshoVortex_3D/run.sh

+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e glassCube_64.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 4 gresho.yml 2>&1 | tee output.log
+
+# Plot the solution
+python plotSolution.py 11

examples/HydrostaticHalo/hydrostatic.yml

 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.9885e39     # 10^6 solar masses
+  UnitMass_in_cgs:     1.98848e39    # 10^6 solar masses
   UnitLength_in_cgs:   3.0856776e21  # Kiloparsecs
   UnitVelocity_in_cgs: 1e5           # Kilometres per second
   UnitCurrent_in_cgs:  1   # Amperes

examples/InteractingBlastWaves_1D/getReference.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/interactingBlastWaves1D_exact.txt