Merge branch 'master' into dopair-vectorisation

Conflicts:
	src/Makefile.am
	src/runner_doiact_vec.c

.gitignore

@@ -76,6 +76,9 @@ tests/testRiemannExact
 tests/testRiemannTRRS
 tests/testRiemannHLLC
 tests/testMatrixInversion
+tests/testDump
+tests/testLogger
+tests/benchmarkInteractions
 
 theory/latex/swift.pdf
 theory/SPH/Kernels/kernels.pdf

configure.ac

@@ -351,7 +351,7 @@ AC_ARG_WITH([tcmalloc],
    [with_tcmalloc="no"]
 )
 if test "x$with_tcmalloc" != "xno"; then
-   if test "x$with_tcmalloc" != "xyes" && test "x$with_tcmalloc" != "x"; then
+   if test "x$with_tcmalloc" != "xyes" -a "x$with_tcmalloc" != "x"; then
       tclibs="-L$with_tcmalloc -ltcmalloc"
    else
       tclibs="-ltcmalloc"
@@ -361,7 +361,7 @@ if test "x$with_tcmalloc" != "xno"; then
 
    #  Could just have the minimal version.
    if test "$have_tcmalloc" = "no"; then
-      if test "x$with_tcmalloc" != "xyes" && test "x$with_tcmalloc" != "x"; then
+      if test "x$with_tcmalloc" != "xyes" -a "x$with_tcmalloc" != "x"; then
          tclibs="-L$with_tcmalloc -ltcmalloc_minimal"
       else
          tclibs="-ltcmalloc_minimal"
@@ -394,7 +394,7 @@ AC_ARG_WITH([profiler],
    [with_profiler="yes"]
 )
 if test "x$with_profiler" != "xno"; then
-   if test "x$with_profiler" != "xyes" && test "x$with_profiler" != "x"; then
+   if test "x$with_profiler" != "xyes" -a "x$with_profiler" != "x"; then
       proflibs="-L$with_profiler -lprofiler"
    else
       proflibs="-lprofiler"
@@ -411,6 +411,38 @@ fi
 AC_SUBST([PROFILER_LIBS])
 AM_CONDITIONAL([HAVEPROFILER],[test -n "$PROFILER_LIBS"])
 
+#  Check for jemalloc another fast malloc that is good with contention.
+have_jemalloc="no"
+AC_ARG_WITH([jemalloc],
+   [AS_HELP_STRING([--with-jemalloc],
+      [use jemalloc library or specify the directory with lib @<:@yes/no@:>@]
+   )],
+   [with_jemalloc="$withval"],
+   [with_jemalloc="no"]
+)
+if test "x$with_jemalloc" != "xno"; then
+   if test "x$with_jemalloc" != "xyes" -a "x$with_jemalloc" != "x"; then
+      jelibs="-L$with_jemalloc -ljemalloc"
+   else
+      jelibs="-ljemalloc"
+   fi
+   AC_CHECK_LIB([jemalloc],[malloc_usable_size],[have_jemalloc="yes"],[have_jemalloc="no"],
+                $jelibs)
+
+   if test "$have_jemalloc" = "yes"; then
+      JEMALLOC_LIBS="$jelibs"
+   else
+      JEMALLOC_LIBS=""
+   fi
+fi
+AC_SUBST([JEMALLOC_LIBS])
+AM_CONDITIONAL([HAVEJEMALLOC],[test -n "$JEMALLOC_LIBS"])
+
+#  Don't allow both tcmalloc and jemalloc.
+if test "x$have_tcmalloc" != "xno" -a "x$have_jemalloc" != "xno"; then
+   AC_MSG_ERROR([Cannot use tcmalloc at same time as jemalloc])
+fi
+
 # Check for HDF5. This is required.
 AX_LIB_HDF5
 
@@ -734,9 +766,6 @@ case "$with_potential" in
    isothermal)
       AC_DEFINE([EXTERNAL_POTENTIAL_ISOTHERMAL], [1], [Isothermal external potential])
    ;; 
-   softened-isothermal)
-      AC_DEFINE([EXTERNAL_POTENTIAL_SOFTENED_ISOTHERMAL], [1], [Softened isothermal external potential])
-   ;; 
    disc-patch)
       AC_DEFINE([EXTERNAL_POTENTIAL_DISC_PATCH], [1], [Disc-patch external potential])
    ;; 
@@ -781,6 +810,7 @@ AC_MSG_RESULT([
    FFTW3 enabled   : $have_fftw3
    libNUMA enabled : $have_numa
    Using tcmalloc  : $have_tcmalloc
+   Using jemalloc  : $have_jemalloc
    CPU profiler    : $have_profiler
 
    Hydro scheme       : $with_hydro
@@ -795,5 +825,8 @@ AC_MSG_RESULT([
    Debugging checks   : $enable_debugging_checks
 ])
 
+# Make sure the latest git revision string gets included
+touch src/version.c
+
 # Generate output.
 AC_OUTPUT

examples/CoolingBox/coolingBox.yml

 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     2.0e33   # Solar masses
-  UnitLength_in_cgs:   3.0857e21   # Kiloparsecs
-  UnitVelocity_in_cgs: 1.0e5   # Time unit is cooling time
-  UnitCurrent_in_cgs:  1   # Amperes
-  UnitTemp_in_cgs:     1   # Kelvin
+  UnitMass_in_cgs:     2.0e33     # Solar masses
+  UnitLength_in_cgs:   3.0857e21  # Kiloparsecs
+  UnitVelocity_in_cgs: 1.0e5      # Kilometers 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:   4.    # The end time of the simulation (in internal units).
-  dt_min:     1e-4  # 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).
+  time_end:   0.25  # The end time of the simulation (in internal units).
+  dt_min:     1e-5  # 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:            coolingBox # Common part of the name of output files
   time_first:          0.         # Time of the first output (in internal units)
-  delta_time:          1.0e-1       # Time difference between consecutive outputs (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-2 # Time between statistics output
+  delta_time:          1e-3 # Time between statistics output
 
 # Parameters for the hydrodynamics scheme
 SPH:
@@ -35,8 +35,8 @@ InitialConditions:
 
 # Dimensionless pre-factor for the time-step condition
 LambdaCooling:
-  lambda_cgs:                 1.0e-22    # Cooling rate (in cgs units)
-  minimum_temperature:         1.0e4  # Minimal temperature (Kelvin)
-  mean_molecular_weight:       0.59   # Mean molecular weight
-  hydrogen_mass_abundance:     0.75   # Hydrogen mass abundance (dimensionless)
-  cooling_tstep_mult:          1.0    # Dimensionless pre-factor for the time-step condition
+  lambda_cgs:                  1.0e-22    # Cooling rate (in cgs units)
+  minimum_temperature:         1.0e4      # Minimal temperature (Kelvin)
+  mean_molecular_weight:       0.59       # Mean molecular weight
+  hydrogen_mass_abundance:     0.75       # Hydrogen mass abundance (dimensionless)
+  cooling_tstep_mult:          1.0        # Dimensionless pre-factor for the time-step condition

examples/CoolingBox/energy_plot.py

+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' : (3.15,3.15),
+'figure.subplot.left'    : 0.145,
+'figure.subplot.right'   : 0.99,
+'figure.subplot.bottom'  : 0.11,
+'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']})
+
+
 import numpy as np
-import matplotlib.pyplot as plt
 import h5py as h5
 import sys
 
+# File containing the total energy
 stats_filename = "./energy.txt"
+
+# First snapshot
 snap_filename = "coolingBox_000.hdf5"
-#plot_dir = "./"
 
-#some constants in cgs units
+# Some constants in cgs units
 k_b = 1.38E-16 #boltzmann
 m_p = 1.67e-24 #proton mass
-#initial conditions set in makeIC.py
-rho = 3.2e3
-P = 4.5e6
-#n_H_cgs = 0.0001
-gamma = 5./3.
+
+# Initial conditions set in makeIC.py
 T_init = 1.0e5
 
-#Read the units parameters from the snapshot
+# Read the initial state of the gas
 f = h5.File(snap_filename,'r')
+rho = np.mean(f["/PartType0/Density"])
+pressure = np.mean(f["/PartType0/Pressure"])
+
+# Read the units parameters from the snapshot
 units = f["InternalCodeUnits"]
 unit_mass = units.attrs["Unit mass in cgs (U_M)"]
 unit_length = units.attrs["Unit length in cgs (U_L)"]
 unit_time = units.attrs["Unit time in cgs (U_t)"]
+
+# Read the properties of the cooling function
 parameters = f["Parameters"]
 cooling_lambda = float(parameters.attrs["LambdaCooling:lambda_cgs"])
 min_T = float(parameters.attrs["LambdaCooling:minimum_temperature"])
 mu = float(parameters.attrs["LambdaCooling:mean_molecular_weight"])
 X_H = float(parameters.attrs["LambdaCooling:hydrogen_mass_abundance"])
 
-#get number of particles
-header = f["Header"]
-n_particles = header.attrs["NumPart_ThisFile"][0]
+# Read the adiabatic index
+gamma = float(f["HydroScheme"].attrs["Adiabatic index"])
+
+print "Initial density :", rho
+print "Initial pressure:", pressure
+print "Adiabatic index :", gamma
 
-#read energy and time arrays
+# Read energy and time arrays
 array = np.genfromtxt(stats_filename,skip_header = 1)
 time = array[:,0]
-kin_plus_therm = array[:,2]
-radiated = array[:,6]
 total_mass = array[:,1]
-
-#ignore first row where there are just zeros
-time = time[1:]
-kin_plus_therm = kin_plus_therm[1:]
-radiated = radiated[1:]
-total_mass = total_mass[1:]
-
-total_energy = kin_plus_therm + radiated
+total_energy = array[:,2]
+kinetic_energy = array[:,3]
+internal_energy = array[:,4]
+radiated_energy = array[:,8]
 initial_energy = total_energy[0]
-#conversions to cgs
+
+# Conversions to cgs
 rho_cgs = rho * unit_mass / (unit_length)**3
 time_cgs = time * unit_time
-initial_energy_cgs = initial_energy/total_mass[0] * unit_length**2 / (unit_time)**2 
-n_H_cgs = X_H * rho_cgs / m_p
+total_energy_cgs = total_energy / total_mass[0] * unit_length**2 / (unit_time)**2
+kinetic_energy_cgs = kinetic_energy / total_mass[0] * unit_length**2 / (unit_time)**2
+internal_energy_cgs = internal_energy / total_mass[0] * unit_length**2 / (unit_time)**2
+radiated_energy_cgs = radiated_energy / total_mass[0] * unit_length**2 / (unit_time)**2  
 
-#find the energy floor
+# Find the energy floor
 u_floor_cgs = k_b * min_T / (mu * m_p * (gamma - 1.))
 
-#find analytic solution
-analytic_time_cgs = np.linspace(0,time_cgs[-1],1000)
+# Find analytic solution
+initial_energy_cgs = initial_energy/total_mass[0] * unit_length**2 / (unit_time)**2 
+n_H_cgs = X_H * rho_cgs / m_p
 du_dt_cgs = -cooling_lambda * n_H_cgs**2 / rho_cgs
+cooling_time_cgs = (initial_energy_cgs/(-du_dt_cgs))[0]
+analytic_time_cgs = np.linspace(0, cooling_time_cgs * 1.8, 1000)
 u_analytic_cgs = du_dt_cgs*analytic_time_cgs + initial_energy_cgs
-cooling_time_cgs = initial_energy_cgs/(-du_dt_cgs)
-
-for i in range(u_analytic_cgs.size):
-    if u_analytic_cgs[i]<u_floor_cgs:
-        u_analytic_cgs[i] = u_floor_cgs
-
-#rescale analytic solution
-u_analytic = u_analytic_cgs/initial_energy_cgs
-
-#put time in units of cooling_time
-time=time_cgs/cooling_time_cgs
-analytic_time = analytic_time_cgs/cooling_time_cgs
-
-#rescale (numerical) energy by initial energy
-radiated /= initial_energy
-kin_plus_therm /= initial_energy
-total_energy = kin_plus_therm + radiated
-plt.plot(time,kin_plus_therm,'kd',label = "Kinetic + thermal energy")
-plt.plot(time,radiated,'bo',label = "Radiated energy")
-plt.plot(time,total_energy,'g',label = "Total energy")
-plt.plot(analytic_time,u_analytic,'r',lw = 2.0,label = "Analytic Solution")
-#plt.plot(analytic_time,1-u_analytic,'k',lw = 2.0)
-#plt.plot((cooling_time,cooling_time),(0,1),'b',label = "Cooling time")
-#plt.plot((time[1]-time_cgs[0],time_cgs[1]-time_cgs[0]),(0,1),'m',label = "First output")
-#plt.title(r"$n_H = %1.1e \, \mathrm{cm}^{-3}$" %n_H_cgs)
-plt.xlabel("Time / cooling time")
-plt.ylabel("Energy / Initial energy")
-#plt.ylim(0,1.1)
-plt.ylim(0.999,1.001)
-#plt.xlim(0,min(10,time[-1]))
-plt.legend(loc = "upper right")    
-if (int(sys.argv[1])==0):
-    plt.show()
-else:
-    plt.savefig(full_plot_filename,format = "png")
-    plt.close()
+u_analytic_cgs[u_analytic_cgs < u_floor_cgs] = u_floor_cgs
+
+print "Cooling time:", cooling_time_cgs, "[s]"
+
+# Read snapshots
+u_snapshots_cgs = zeros(25)
+t_snapshots_cgs = zeros(25)
+for i in range(25):
+    snap = h5.File("coolingBox_%0.3d.hdf5"%i,'r')
+    u_snapshots_cgs[i] = sum(snap["/PartType0/InternalEnergy"][:] * snap["/PartType0/Masses"][:])  / total_mass[0] * unit_length**2 / (unit_time)**2
+    t_snapshots_cgs[i] = snap["/Header"].attrs["Time"] * unit_time
+
+
+figure()
+plot(time_cgs, total_energy_cgs, 'r-', lw=1.6, label="Gas total energy")
+plot(t_snapshots_cgs, u_snapshots_cgs, 'rD', ms=3)
+plot(time_cgs, radiated_energy_cgs, 'g-', lw=1.6, label="Radiated energy")
+plot(time_cgs, total_energy_cgs + radiated_energy_cgs, 'b-', lw=0.6, label="Gas total + radiated")
+
+plot(analytic_time_cgs, u_analytic_cgs, '--', color='k', alpha=0.8, lw=1.0, label="Analytic solution")
+
+legend(loc="upper right", fontsize=8, frameon=False, handlelength=3, ncol=1)
+xlabel("${\\rm{Time~[s]}}$", labelpad=0)
+ylabel("${\\rm{Energy~[erg]}}$")
+xlim(0, 1.5*cooling_time_cgs)
+ylim(0, 1.5*u_analytic_cgs[0])
+
+savefig("energy.png", dpi=200)
+
+

examples/CoolingBox/getGlass.sh

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

examples/CoolingBox/makeIC.py

 ###############################################################################
  # This file is part of SWIFT.
- # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@durhama.ac.uk)
  #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  # 
  # This program is free software: you can redistribute it and/or modify
@@ -22,13 +22,11 @@ import h5py
 import sys
 from numpy import *
 
-# Generates a swift IC file containing a cartesian distribution of particles
-# at a constant density and pressure in a cubic box
+# Generates a SWIFT IC file with a constant density and pressure
 
 # Parameters
 periodic= 1           # 1 For periodic box
 boxSize = 1           # 1 kiloparsec    
-L = int(sys.argv[1])  # Number of particles along one axis
 rho = 3.2e3           # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
 P = 4.5e6             # Pressure in code units (at 10^5K)
 gamma = 5./3.         # Gas adiabatic index
@@ -36,12 +34,17 @@ eta = 1.2349          # 48 ngbs with cubic spline kernel
 fileName = "coolingBox.hdf5" 
 
 #---------------------------------------------------
-numPart = L**3
-mass = boxSize**3 * rho / numPart
-print mass
-internalEnergy = P / ((gamma - 1.)*rho)
 
-#--------------------------------------------------
+# Read id, position and h from glass
+glass = h5py.File("glassCube_32.hdf5", "r")
+ids = glass["/PartType0/ParticleIDs"][:]
+pos = glass["/PartType0/Coordinates"][:,:] * boxSize
+h = glass["/PartType0/SmoothingLength"][:] * boxSize
+
+# Compute basic properties
+numPart = size(pos) / 3
+mass = boxSize**3 * rho / numPart
+internalEnergy = P / ((gamma - 1.) * rho)
 
 #File
 file = h5py.File(fileName, 'w')
@@ -57,11 +60,11 @@ 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
+# Runtime parameters
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
-#Units
+# Units
 grp = file.create_group("/Units")
 grp.attrs["Unit length in cgs (U_L)"] = 3.0857e21 
 grp.attrs["Unit mass in cgs (U_M)"] = 2.0e33 
@@ -75,35 +78,26 @@ grp = file.create_group("/PartType0")
 v  = zeros((numPart, 3))
 ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
 ds[()] = v
-v = zeros(1)
 
 m = full((numPart, 1), mass)
 ds = grp.create_dataset('Masses', (numPart,1), 'f')
 ds[()] = m
-m = zeros(1)
 
-h = full((numPart, 1), eta * boxSize / L)
-ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
+h = reshape(h, (numPart, 1))
+ds = grp.create_dataset('SmoothingLength', (numPart, 1), 'f')
 ds[()] = h
-h = zeros(1)
 
 u = full((numPart, 1), internalEnergy)
 ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
 ds[()] = u
-u = zeros(1)
 
-
-ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
+ids = reshape(ids, (numPart, 1))
 ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
-ds[()] = ids + 1
-x      = ids % L;
-y      = ((ids - x) / L) % L;
-z      = (ids - x - L * y) / L**2;
-coords = zeros((numPart, 3))
-coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
-coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
-coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
+ds[()] = ids
+
 ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
-ds[()] = coords
+ds[()] = pos
 
 file.close()
+
+print numPart

examples/CoolingBox/run.sh

+
 #!/bin/bash
 
 # Generate the initial conditions if they are not present.
-echo "Generating initial conditions for the cooling box example..."
-
-python makeIC.py 10
-
-../swift -s -C -t 16 coolingBox.yml 
-
-#-C 2>&1 | tee output.log
+if [ ! -e glassCube_32.hdf5 ]
+then
+    echo "Fetching initial glass file for the cooling box example..."
+    ./getGlass.sh
+fi
+if [ ! -e coolingBox.hdf5 ]
+then
+    echo "Generating initial conditions for the cooling box example..."
+    python makeIC.py
+fi
 
-python energy_plot.py 0
+# Run SWIFT
+../swift -s -C -t 1 coolingBox.yml
 
-#python test_energy_conservation.py 0
+# Check energy conservation and cooling rate
+python energy_plot.py

examples/CoolingBox/test_energy_conservation.py

-import numpy as np
-import matplotlib.pyplot as plt
-import h5py as h5
-import sys
-
-stats_filename = "./energy.txt"
-snap_filename = "coolingBox_000.hdf5"
-#plot_dir = "./"
-n_snaps = 41
-time_end = 4.0
-dt_snap = 0.1
-#some constants in cgs units
-k_b = 1.38E-16 #boltzmann
-m_p = 1.67e-24 #proton mass
-#initial conditions set in makeIC.py
-rho = 4.8e3
-P = 4.5e6
-#n_H_cgs = 0.0001
-gamma = 5./3.
-T_init = 1.0e5
-
-#find the sound speed
-
-#Read the units parameters from the snapshot
-f = h5.File(snap_filename,'r')
-units = f["InternalCodeUnits"]
-unit_mass = units.attrs["Unit mass in cgs (U_M)"]
-unit_length = units.attrs["Unit length in cgs (U_L)"]
-unit_time = units.attrs["Unit time in cgs (U_t)"]
-parameters = f["Parameters"]
-cooling_lambda = float(parameters.attrs["LambdaCooling:lambda_cgs"])
-min_T = float(parameters.attrs["LambdaCooling:minimum_temperature"])
-mu = float(parameters.attrs["LambdaCooling:mean_molecular_weight"])
-X_H = float(parameters.attrs["LambdaCooling:hydrogen_mass_abundance"])
-
-#get number of particles
-header = f["Header"]
-n_particles = header.attrs["NumPart_ThisFile"][0]
-#read energy and time arrays
-array = np.genfromtxt(stats_filename,skip_header = 1)
-time = array[:,0]
-total_energy = array[:,2]
-total_mass = array[:,1]
-
-time = time[1:]
-total_energy = total_energy[1:]
-total_mass = total_mass[1:]
-
-#conversions to cgs
-rho_cgs = rho * unit_mass / (unit_length)**3
-time_cgs = time * unit_time
-u_init_cgs = total_energy[0]/(total_mass[0]) * unit_length**2 / (unit_time)**2 
-n_H_cgs = X_H * rho_cgs / m_p
-
-#find the sound speed in cgs
-c_s = np.sqrt((gamma - 1.)*k_b*T_init/(mu*m_p))
-#assume box size is unit length
-sound_crossing_time = unit_length/c_s
-
-print "Sound speed = %g cm/s" %c_s
-print "Sound crossing time = %g s" %sound_crossing_time 
-#find the energy floor
-u_floor_cgs = k_b * min_T / (mu * m_p * (gamma - 1.))
-#find analytic solution
-analytic_time_cgs = np.linspace(time_cgs[0],time_cgs[-1],1000)
-du_dt_cgs = -cooling_lambda * n_H_cgs**2 / rho_cgs
-u_analytic = du_dt_cgs*(analytic_time_cgs - analytic_time_cgs[0]) + u_init_cgs
-cooling_time = u_init_cgs/(-du_dt_cgs)
-
-#put time in units of sound crossing time
-time=time_cgs/sound_crossing_time
-analytic_time = analytic_time_cgs/sound_crossing_time
-#rescale energy to initial energy
-total_energy /= total_energy[0]
-u_analytic /= u_init_cgs
-u_floor_cgs /= u_init_cgs
-# plot_title = r"$\Lambda \, = \, %1.1g \mathrm{erg}\mathrm{cm^3}\mathrm{s^{-1}} \, \, T_{init} = %1.1g\mathrm{K} \, \, T_{floor} = %1.1g\mathrm{K} \, \, n_H = %1.1g\mathrm{cm^{-3}}$" %(cooling_lambda,T_init,T_floor,n_H)
-# plot_filename = "energy_plot_creasey_no_cooling_T_init_1p0e5_n_H_0p1.png"
-#analytic_solution = np.zeros(n_snaps-1)
-for i in range(u_analytic.size):
-    if u_analytic[i]<u_floor_cgs:
-        u_analytic[i] = u_floor_cgs
-plt.plot(time-time[0],total_energy,'k',label = "Numerical solution from energy.txt")
-plt.plot(analytic_time-analytic_time[0],u_analytic,'r',lw = 2.0,label = "Analytic Solution")
-
-#now get energies from the snapshots
-snapshot_time = np.linspace(0,time_end,num = n_snaps)
-snapshot_time = snapshot_time[1:]
-snapshot_time_cgs = snapshot_time * unit_time
-snapshot_time = snapshot_time_cgs/ sound_crossing_time
-snapshot_time -= snapshot_time[0]
-snapshot_energy = np.zeros(n_snaps)
-for i in range(0,n_snaps):
-    snap_filename = "coolingBox_%03d.hdf5" %i
-    f = h5.File(snap_filename,'r')
-    snapshot_internal_energy_array = np.array(f["PartType0/InternalEnergy"])
-    total_internal_energy = np.sum(snapshot_internal_energy_array)
-    velocity_array = np.array(f["PartType0/Velocities"])
-    total_kinetic_energy = 0.5*np.sum(velocity_array**2)
-    snapshot_energy[i] = total_internal_energy + total_kinetic_energy
-snapshot_energy/=snapshot_energy[0]
-snapshot_energy = snapshot_energy[1:]
-
-plt.plot(snapshot_time,snapshot_energy,'bd',label = "Numerical solution from snapshots")
-
-#plt.title(r"$n_H = %1.1e \, \mathrm{cm}^{-3}$" %n_H_cgs)
-plt.xlabel("Time (sound crossing time)")
-plt.ylabel("Energy/Initial energy")
-plt.ylim(0.99,1.01)
-#plt.xlim(0,min(10,time[-1]))
-plt.legend(loc = "upper right")    
-if (int(sys.argv[1])==0):
-    plt.show()
-else:
-    plt.savefig(full_plot_filename,format = "png")
-    plt.close()

examples/CoolingHalo/cooling_halo.yml

@@ -37,7 +37,7 @@ InitialConditions:
   shift_z:    0.
  
 # External potential parameters
-SoftenedIsothermalPotential:
+IsothermalPotential:
   position_x:      0.     # location of centre of isothermal potential in internal units
   position_y:      0.
   position_z:      0.	

examples/CoolingHaloWithSpin/cooling_halo.yml

@@ -9,8 +9,8 @@ InternalUnitSystem:
 # 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-7  # The minimal time-step size of the simulation (in internal units).
+  time_end:   10.   # The end time of the simulation (in internal units).
+  dt_min:     1e-5  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-1  # The maximal time-step size of the simulation (in internal units).
 
 # Parameters governing the conserved quantities statistics
@@ -34,13 +34,13 @@ InitialConditions:
   file_name:  CoolingHalo.hdf5       # The file to read
  
 # External potential parameters
-SoftenedIsothermalPotential:
-  position_x:      0.     # location of centre of isothermal potential in internal units
+IsothermalPotential:
+  position_x:      0.     # Location of centre of isothermal potential in internal units
   position_y:      0.
   position_z:      0.	
-  vrot:            200.     # rotation speed of isothermal potential in internal units
-  timestep_mult:   0.03     # controls time step
-  epsilon:         1.0      #softening for the isothermal potential
+  vrot:            200.   # Rotation speed of isothermal potential in internal units
+  timestep_mult:   0.03   # Controls time step
+  epsilon:         1.0    # Softening for the isothermal potential
 
 # Cooling parameters
 LambdaCooling:

examples/CoolingHaloWithSpin/density_profile.py

@@ -28,7 +28,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
 lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
 X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
@@ -101,18 +101,18 @@ for i in range(n_snaps):
         rho_0 = density[0]
 
         rho_analytic_init = rho_0 * (radial_bin_mids/r_0)**(-2)
-    plt.plot(radial_bin_mids,density/rho_analytic_init,'ko',label = "Average density of shell")
-    #plt.plot(t,rho_analytic,label = "Initial analytic density profile"
+    plt.plot(radial_bin_mids,density,'ko',label = "Average density of shell")
+    plt.plot(radial_bin_mids,rho_analytic_init,label = "Initial analytic density profile")
     plt.xlabel(r"$r / r_{vir}$")
-    plt.ylabel(r"$\rho / \rho_{init})$")
+    plt.ylabel(r"$(\rho / \rho_{init})$")
     plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
     #plt.ylim((1.e-2,1.e1))
-    plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,20),'r',label = "Cooling radius")
+    plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(1.0e-4,1.0e4),'r',label = "Cooling radius")
     plt.xlim((radial_bin_mids[0],max_r))
-    plt.ylim((0,20))
-    plt.plot((0,max_r),(1,1))
+    plt.ylim((1.0e-4,1.0e4))
+    #plt.plot((0,max_r),(1,1))
     #plt.xscale('log')
-    #plt.yscale('log')
+    plt.yscale('log')
     plt.legend(loc = "upper right")
     plot_filename = "density_profile_%03d.png" %i
     plt.savefig(plot_filename,format = "png")

examples/CoolingHaloWithSpin/internal_energy_profile.py

@@ -46,7 +46,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
 lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
 X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])

examples/CoolingHaloWithSpin/makeIC.py

@@ -36,6 +36,7 @@ h = 0.67777 # hubble parameter
 gamma = 5./3.
 eta = 1.2349
 spin_lambda = 0.05 #spin parameter
+f_b = 0.2 #baryon fraction
 
 # First set unit velocity and then the circular velocity parameter for the isothermal potential
 const_unit_velocity_in_cgs = 1.e5 #kms^-1
@@ -99,6 +100,8 @@ grp.attrs["PeriodicBoundariesOn"] = periodic
 # set seed for random number
 np.random.seed(1234)
 
+gas_mass = f_b * np.sqrt(3.) / 2. #virial mass of halo is 1, virial radius is 1, enclosed mass scales with r
+gas_particle_mass = gas_mass / float(N)
 
 # Positions
 # r^(-2) distribution corresponds to uniform distribution in radius
@@ -164,12 +167,12 @@ N = x_coords.size
 print "Number of particles in the box = " , N
 
 #make the coords and radius arrays again
-coords_2 = np.zeros((N,3))
-coords_2[:,0] = x_coords
-coords_2[:,1] = y_coords
-coords_2[:,2] = z_coords
+coords= np.zeros((N,3))
+coords[:,0] = x_coords
+coords[:,1] = y_coords
+coords[:,2] = z_coords
 
-radius = np.sqrt((coords_2[:,0]-boxSize/2.)**2 + (coords_2[:,1]-boxSize/2.)**2 + (coords_2[:,2]-boxSize/2.)**2)
+radius = np.sqrt((coords[:,0]-boxSize/2.)**2 + (coords[:,1]-boxSize/2.)**2 + (coords[:,2]-boxSize/2.)**2)
 
 #now give particle's velocities
 v = np.zeros((N,3))
@@ -184,7 +187,7 @@ print "J =", J
 omega = np.zeros((N,3))
 for i in range(N):
     omega[i,2] = 3.*J / radius[i]
-    v[i,:] = np.cross(omega[i,:],(coords_2[i,:]-boxSize/2.))
+    v[i,:] = np.cross(omega[i,:],(coords[i,:]-boxSize/2.))
         
 # Header
 grp = file.create_group("/Header")
@@ -202,16 +205,15 @@ grp.attrs["Dimension"] = 3
 grp = file.create_group("/PartType0")
 
 ds = grp.create_dataset('Coordinates', (N, 3), 'd')
-ds[()] = coords_2
-coords_2 = np.zeros(1)
+ds[()] = coords
+coords = np.zeros(1)
 
 ds = grp.create_dataset('Velocities', (N, 3), 'f')
 ds[()] = v
 v = np.zeros(1)
 
 # All particles of equal mass
-mass = 1. / N
-m = np.full((N,),mass)
+m = np.full((N,),gas_particle_mass)
 ds = grp.create_dataset('Masses', (N, ), 'f')
 ds[()] = m
 m = np.zeros(1)

examples/CoolingHaloWithSpin/run.sh

@@ -4,7 +4,8 @@
 echo "Generating initial conditions for the isothermal potential box example..."
 python makeIC.py 10000 
 
-../swift -g -s -C -t 16 cooling_halo.yml 2>&1 | tee output.log
+# Run SWIFT with external potential, SPH and cooling
+../swift -g -s -C -t 1 cooling_halo.yml 2>&1 | tee output.log
 
 # python radial_profile.py 10
 

examples/CoolingHaloWithSpin/velocity_profile.py

@@ -46,7 +46,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
 header = f["Header"]
 N = header.attrs["NumPart_Total"][0]

examples/EAGLE_100/README

+ICs extracted from the EAGLE suite of simulations. 
+
+WARNING: The ICs are 217GB in size. They contain ~3.4G DM particles,
+~3.2G gas particles and ~170M star particles
+
+The particle distribution here is the snapshot 27 (z=0.1) of the 100Mpc
+Ref-model. h- and a- factors from the original Gadget code have been
+corrected for. Variables not used in a pure hydro & gravity code have
+been removed. 
+Everything is ready to be run without cosmological integration. 
+
+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

examples/Feedback/feedback.yml → examples/EAGLE_100/eagle_100.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
+  UnitMass_in_cgs:     1.989e43      # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in 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).
+  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_max:     1e-4  # The maximal time-step size of the simulation (in internal units).
 
 # Parameters governing the snapshots
 Snapshots:
-  basename:            Feedback # 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)
+  basename:            eagle # Common part of the name of output files
+  time_first:          0.    # Time of the first output (in internal units)
+  delta_time:          1e-3  # Time difference between consecutive outputs (in internal units)
 
 # Parameters governing the conserved quantities statistics
 Statistics:
-  delta_time:          1e-3 # Time between statistics output
+  delta_time:          1e-2 # Time between statistics output
 
 # Parameters for the hydrodynamics scheme
 SPH:
@@ -31,13 +31,5 @@ SPH:
 
 # Parameters related to the initial conditions
 InitialConditions:
-  file_name:  ./Feedback.hdf5          # The file to read
+  file_name:  ./EAGLE_ICs_100.hdf5     # The file to read
 
-# Parameters for feedback
-
-SN:
-  time:    0.001 # time the SN explodes (internal units)
-  energy:  1.0   # energy of the explosion (internal units)
-  x:       0.5   # x-position of explostion (internal units)
-  y:       0.5   # y-position of explostion (internal units)
-  z:       0.5   # z-position of explostion (internal units)

examples/EAGLE_100/getIC.sh

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

examples/EAGLE_100/run.sh

+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e EAGLE_ICs_100.hdf5 ]
+then
+    echo "Fetching initial conditions for the EAGLE 100Mpc example..."
+    ./getIC.sh
+fi
+
+../swift -s -t 16 eagle_100.yml 2>&1 | tee output.log
+

examples/ExternalPointMass/energy_plot.py

+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' : (3.15,3.15),
+'figure.subplot.left'    : 0.145,
+'figure.subplot.right'   : 0.99,
+'figure.subplot.bottom'  : 0.11,
+'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']})
+
+
+import numpy as np
+import h5py as h5
+import sys
+
+# File containing the total energy
+stats_filename = "./energy.txt"
+
+# First snapshot
+snap_filename = "pointMass_000.hdf5"
+f = h5.File(snap_filename,'r')
+
+# Read the units parameters from the snapshot
+units = f["InternalCodeUnits"]
+unit_mass = units.attrs["Unit mass in cgs (U_M)"]
+unit_length = units.attrs["Unit length in cgs (U_L)"]
+unit_time = units.attrs["Unit time in cgs (U_t)"]
+
+G = 6.67408e-8 * unit_mass * unit_time**2 / unit_length**3
+
+# Read the header
+header = f["Header"]
+box_size = float(header.attrs["BoxSize"][0])
+
+# Read the properties of the potential
+parameters = f["Parameters"]
+mass = float(parameters.attrs["PointMassPotential:mass"])
+centre = [box_size/2, box_size/2, box_size/2]
+f.close()
+
+# Read the statistics summary
+file_energy = np.loadtxt("energy.txt")
+time_stats = file_energy[:,0]
+E_kin_stats = file_energy[:,3]
+E_pot_stats = file_energy[:,5]
+E_tot_stats = E_kin_stats + E_pot_stats
+
+# Read the snapshots
+time_snap = np.zeros(402)
+E_kin_snap = np.zeros(402)
+E_pot_snap = np.zeros(402)
+E_tot_snap = np.zeros(402)
+Lz_snap = np.zeros(402)
+
+# Read all the particles from the snapshots
+for i in range(402):
+    snap_filename = "pointMass_%0.3d.hdf5"%i
+    f = h5.File(snap_filename,'r')
+
+    pos_x = f["PartType1/Coordinates"][:,0]
+    pos_y = f["PartType1/Coordinates"][:,1]
+    pos_z = f["PartType1/Coordinates"][:,2]
+    vel_x = f["PartType1/Velocities"][:,0]
+    vel_y = f["PartType1/Velocities"][:,1]
+    vel_z = f["PartType1/Velocities"][:,2]
+    m = f["/PartType1/Masses"][:]
+    
+    r = np.sqrt((pos_x[:] - centre[0])**2 + (pos_y[:] - centre[1])**2 + (pos_z[:] - centre[2])**2)
+    Lz = (pos_x[:] - centre[0]) * vel_y[:] - (pos_y[:] - centre[1]) * vel_x[:]
+
+    time_snap[i] = f["Header"].attrs["Time"]
+    E_kin_snap[i] = np.sum(0.5 * m * (vel_x[:]**2 + vel_y[:]**2 + vel_z[:]**2))
+    E_pot_snap[i] = np.sum(-mass * m * G / r)
+    E_tot_snap[i] = E_kin_snap[i] + E_pot_snap[i]
+    Lz_snap[i] = np.sum(Lz)
+
+# Plot energy evolution
+figure()
+plot(time_stats, E_kin_stats, "r-", lw=0.5, label="Kinetic energy")
+plot(time_stats, E_pot_stats, "g-", lw=0.5, label="Potential energy")
+plot(time_stats, E_tot_stats, "k-", lw=0.5, label="Total energy")
+
+plot(time_snap[::10], E_kin_snap[::10], "rD", lw=0.5, ms=2)
+plot(time_snap[::10], E_pot_snap[::10], "gD", lw=0.5, ms=2)
+plot(time_snap[::10], E_tot_snap[::10], "kD", lw=0.5, ms=2)
+
+legend(loc="center right", fontsize=8, frameon=False, handlelength=3, ncol=1)
+xlabel("${\\rm{Time}}$", labelpad=0)
+ylabel("${\\rm{Energy}}$",labelpad=0)
+xlim(0, 8)
+
+savefig("energy.png", dpi=200)
+
+# Plot angular momentum evolution
+figure()
+plot(time_snap, Lz_snap, "k-", lw=0.5, ms=2)
+
+xlabel("${\\rm{Time}}$", labelpad=0)
+ylabel("${\\rm{Angular~momentum}}$",labelpad=0)
+xlim(0, 8)
+
+savefig("angular_momentum.png", dpi=200)
+
+