Merge branch 'master' into cbrt

.gitignore

@@ -22,8 +22,6 @@ doc/Doxyfile
 
 examples/swift
 examples/swift_mpi
-examples/swift_fixdt
-examples/swift_fixdt_mpi
 examples/*.xmf
 examples/used_parameters.yml
 examples/energy.txt

README

@@ -13,8 +13,6 @@ See INSTALL.swift for install instructions.
 
 Usage: swift [OPTION]... PARAMFILE
        swift_mpi [OPTION]... PARAMFILE
-       swift_fixdt [OPTION]... PARAMFILE
-       swift_fixdt_mpi [OPTION]... PARAMFILE
 
 Valid options are:
   -a          Pin runners using processor affinity

configure.ac

@@ -170,6 +170,18 @@ if test "x$enable_debug" = "xyes"; then
    fi
 fi
 
+# Check if task debugging is on.
+AC_ARG_ENABLE([task-debugging],
+   [AS_HELP_STRING([--enable-task-debugging],
+     [Store task timing information and generate task dump files @<:@yes/no@:>@]
+   )],
+   [enable_task_debugging="$enableval"],
+   [enable_task_debugging="no"]
+)
+if test "$enable_task_debugging" = "yes"; then
+   AC_DEFINE([SWIFT_DEBUG_TASKS],1,[Enable task debugging])
+fi
+
 # Define HAVE_POSIX_MEMALIGN if it works.
 AX_FUNC_POSIX_MEMALIGN
 
@@ -533,6 +545,7 @@ AC_MSG_RESULT([
    libNUMA enabled : $have_numa
    Using tcmalloc  : $have_tcmalloc
    CPU profiler    : $have_profiler
+   Task debugging  : $enable_task_debugging
 ])
 
 # Generate 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.01e21   # Kilparsecs
+  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
@@ -9,9 +9,9 @@ InternalUnitSystem:
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   1.0    # 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).
+  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).
 
 # Parameters governing the snapshots
 Snapshots:
@@ -27,7 +27,6 @@ Statistics:
 SPH:
   resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
   delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
-  max_smoothing_length:  0.1      # Maximal smoothing length allowed (in internal units).
   CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
   
 # Parameters related to the initial conditions
@@ -36,9 +35,8 @@ InitialConditions:
 
 # Dimensionless pre-factor for the time-step condition
 LambdaCooling:
-  lambda:                      0.0    # Cooling rate (in cgs units)
+  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

@@ -13,7 +13,7 @@ m_p = 1.67e-24 #proton mass
 #initial conditions set in makeIC.py
 rho = 3.2e3
 P = 4.5e6
-n_H_cgs = 0.0001
+#n_H_cgs = 0.0001
 gamma = 5./3.
 T_init = 1.0e5
 
@@ -24,7 +24,7 @@ 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"])
+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"])
@@ -32,50 +32,65 @@ 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]
+kin_plus_therm = array[:,2]
+radiated = array[:,6]
 total_mass = array[:,1]
 
+#ignore first row where there are just zeros
 time = time[1:]
-total_energy = total_energy[1:]
+kin_plus_therm = kin_plus_therm[1:]
+radiated = radiated[1:]
 total_mass = total_mass[1:]
 
+total_energy = kin_plus_therm + radiated
+initial_energy = total_energy[0]
 #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 
+initial_energy_cgs = initial_energy/total_mass[0] * unit_length**2 / (unit_time)**2 
+n_H_cgs = X_H * rho_cgs / m_p
 
 #find the energy floor
-print min_T
 u_floor_cgs = k_b * min_T / (mu * m_p * (gamma - 1.))
+
 #find analytic solution
-analytic_time = np.linspace(time_cgs[0],time_cgs[-1],1000)
-print time_cgs[1]
-print analytic_time[1]
+analytic_time_cgs = np.linspace(0,time_cgs[-1],1000)
 du_dt_cgs = -cooling_lambda * n_H_cgs**2 / rho_cgs
-u_analytic = du_dt_cgs*(analytic_time - analytic_time[0]) + u_init_cgs
-cooling_time = u_init_cgs/(-du_dt_cgs)
-#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_cgs,total_energy,'k',label = "Numerical solution")
-plt.plot(analytic_time,u_analytic,'--r',lw = 2.0,label = "Analytic Solution")
-plt.plot((cooling_time,cooling_time),(0,1),'b',label = "Cooling time")
-plt.plot((time_cgs[0],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 (seconds)")
-plt.ylabel("Energy/Initial energy")
+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.0*cooling_time,time_cgs[-1]))
+#plt.xlim(0,min(10,time[-1]))
 plt.legend(loc = "upper right")    
 if (int(sys.argv[1])==0):
     plt.show()

examples/CoolingBox/makeIC.py

@@ -27,10 +27,10 @@ from numpy import *
 
 # Parameters
 periodic= 1           # 1 For periodic box
-boxSize = 1           #1 kiloparsec    
+boxSize = 1           # 1 kiloparsec    
 L = int(sys.argv[1])  # Number of particles along one axis
-rho = 3.2e3          # Density in code units (0.01 hydrogen atoms per cm^3)
-P = 4.5e6          # Pressure in code units (at 10^5K)
+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
 eta = 1.2349          # 48 ngbs with cubic spline kernel
 fileName = "coolingBox.hdf5" 
@@ -63,9 +63,9 @@ grp.attrs["PeriodicBoundariesOn"] = periodic
 
 #Units
 grp = file.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 3.08e21 
+grp.attrs["Unit length in cgs (U_L)"] = 3.0857e21 
 grp.attrs["Unit mass in cgs (U_M)"] = 2.0e33 
-grp.attrs["Unit time in cgs (U_t)"] = 3.08e16 
+grp.attrs["Unit time in cgs (U_t)"] = 3.0857e16 
 grp.attrs["Unit current in cgs (U_I)"] = 1.
 grp.attrs["Unit temperature in cgs (U_T)"] = 1.
 

examples/CoolingBox/run.sh

@@ -5,6 +5,10 @@ echo "Generating initial conditions for the cooling box example..."
 
 python makeIC.py 10
 
-../swift -s -t 1 coolingBox.yml -C 2>&1 | tee output.log
+../swift -s -C -t 16 coolingBox.yml 
+
+#-C 2>&1 | tee output.log
 
 python energy_plot.py 0
+
+#python test_energy_conservation.py 0

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/README

+
+To make the initial conditions we distribute gas particles randomly in
+a cube with a side length twice that of the virial radius. The density
+profile of the gas is proportional to r^(-2) where r is the distance
+from the centre of the cube.
+
+The parameter v_rot (in makeIC.py and cooling.yml) sets the circular
+velocity of the halo, and by extension, the viral radius, viral mass,
+and the internal energy of the gas such that hydrostatic equilibrium
+is achieved.
+
+While the system is initially in hydrostatic equilibrium, the cooling
+of the gas means that the halo will collapse.
+
+To run this example, make such that the code is compiled with either
+the isothermal potential or softened isothermal potential, and
+'const_lambda' cooling, set in src/const.h. In the latter case, a
+(small) value of epsilon needs to be set in cooling.yml.  0.1 kpc
+should work well.
+
+The plotting scripts produce a plot of the density, internal energy
+and radial velocity profile for each
+snapshot. test_energy_conservation.py shows the evolution of energy
+with time. These can be used to check if the example has run properly.
+

examples/CoolingHalo/cooling_halo.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.9885e39     # 10^6 solar masses
+  UnitLength_in_cgs:   3.0856776e21  # Kiloparsecs
+  UnitVelocity_in_cgs: 1e5           # Kilometres per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   10.0    # 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 conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            CoolingHalo  # 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 for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      1.       # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  CoolingHalo.hdf5       # The file to read
+  shift_x:    0.                  # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    0.
+  shift_z:    0.
+ 
+# External potential parameters
+SoftenedIsothermalPotential:
+  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:         0.1    #softening for the isothermal potential
+
+# Cooling parameters
+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

examples/CoolingHalo/density_profile.py

+import numpy as np
+import h5py as h5
+import matplotlib.pyplot as plt
+import sys
+
+n_snaps = 11
+
+#for the plotting
+#n_radial_bins = int(sys.argv[1])
+
+#some constants
+OMEGA = 0.3 # Cosmological matter fraction at z = 0
+PARSEC_IN_CGS = 3.0856776e18
+KM_PER_SEC_IN_CGS = 1.0e5
+CONST_G_CGS = 6.672e-8
+h = 0.67777 # hubble parameter
+gamma = 5./3.
+eta = 1.2349
+H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+
+#read some header/parameter information from the first snapshot
+
+filename = "Hydrostatic_000.hdf5"
+f = h5.File(filename,'r')
+params = f["Parameters"]
+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["IsothermalPotential:vrot"])
+v_c_cgs = v_c * unit_velocity_cgs
+header = f["Header"]
+N = header.attrs["NumPart_Total"][0]
+box_centre = np.array(header.attrs["BoxSize"])
+
+#calculate r_vir and M_vir from v_c
+r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
+M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+
+for i in range(n_snaps):
+
+    filename = "Hydrostatic_%03d.hdf5" %i
+    f = h5.File(filename,'r')
+    coords_dset = f["PartType0/Coordinates"]
+    coords = np.array(coords_dset)
+#translate coords by centre of box
+    header = f["Header"]
+    snap_time = header.attrs["Time"]
+    snap_time_cgs = snap_time * unit_time_cgs
+    coords[:,0] -= box_centre[0]/2.
+    coords[:,1] -= box_centre[1]/2.
+    coords[:,2] -= box_centre[2]/2.
+    radius = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2)
+    radius_cgs = radius*unit_length_cgs
+    radius_over_virial_radius = radius_cgs / r_vir_cgs
+
+    r = radius_over_virial_radius
+
+    # bin_width = 1./n_radial_bins
+#     hist = np.histogram(r,bins = n_radial_bins)[0] # number of particles in each bin
+
+# #find the mass in each radial bin
+
+#     mass_dset = f["PartType0/Masses"]
+# #mass of each particles should be equal
+#     part_mass = np.array(mass_dset)[0]
+#     part_mass_cgs = part_mass * unit_mass_cgs
+#     part_mass_over_virial_mass = part_mass_cgs / M_vir_cgs 
+
+#     mass_hist = hist * part_mass_over_virial_mass
+#     radial_bin_mids = np.linspace(bin_width/2.,1 - bin_width/2.,n_radial_bins)
+# #volume in each radial bin
+#     volume = 4.*np.pi * radial_bin_mids**2 * bin_width
+
+# #now divide hist by the volume so we have a density in each bin
+
+#     density = mass_hist / volume
+
+    # read the densities
+
+    density_dset = f["PartType0/Density"]
+    density = np.array(density_dset)
+    density_cgs = density * unit_mass_cgs / unit_length_cgs**3
+    rho = density_cgs * r_vir_cgs**3 / M_vir_cgs
+
+    t = np.linspace(0.01,2.0,1000)
+    rho_analytic = t**(-2)/(4.*np.pi)
+
+    plt.plot(r,rho,'x',label = "Numerical solution")
+    plt.plot(t,rho_analytic,label = "Analytic Solution")
+    plt.legend(loc = "upper right")
+    plt.xlabel(r"$r / r_{vir}$")
+    plt.ylabel(r"$\rho / (M_{vir} / r_{vir}^3)$")
+    plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
+    #plt.ylim((0.1,40))
+    plt.xscale('log')
+    plt.yscale('log')
+    plot_filename = "density_profile_%03d.png" %i
+    plt.savefig(plot_filename,format = "png")
+    plt.close()
+

examples/CoolingHalo/internal_energy_profile.py

+import numpy as np
+import h5py as h5
+import matplotlib.pyplot as plt
+import sys
+
+def do_binning(x,y,x_bin_edges):
+
+    #x and y are arrays, where y = f(x)
+    #returns number of elements of x in each bin, and the total of the y elements corresponding to those x values
+
+    n_bins = x_bin_edges.size - 1
+    count = np.zeros(n_bins)
+    y_totals = np.zeros(n_bins)
+    
+    for i in range(n_bins):
+        ind = np.intersect1d(np.where(x > bin_edges[i])[0],np.where(x <= bin_edges[i+1])[0])
+        count[i] = ind.size
+        binned_y = y[ind]
+        y_totals[i] = np.sum(binned_y)
+
+    return(count,y_totals)
+
+
+n_snaps = 100
+
+#for the plotting
+n_radial_bins = int(sys.argv[1])
+
+#some constants
+OMEGA = 0.3 # Cosmological matter fraction at z = 0
+PARSEC_IN_CGS = 3.0856776e18
+KM_PER_SEC_IN_CGS = 1.0e5
+CONST_G_CGS = 6.672e-8
+h = 0.67777 # hubble parameter
+gamma = 5./3.
+eta = 1.2349
+H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+
+#read some header/parameter information from the first snapshot
+
+filename = "Hydrostatic_000.hdf5"
+f = h5.File(filename,'r')
+params = f["Parameters"]
+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["IsothermalPotential:vrot"])
+v_c_cgs = v_c * unit_velocity_cgs
+header = f["Header"]
+N = header.attrs["NumPart_Total"][0]
+box_centre = np.array(header.attrs["BoxSize"])
+
+#calculate r_vir and M_vir from v_c
+r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
+M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+
+for i in range(n_snaps):
+
+    filename = "Hydrostatic_%03d.hdf5" %i
+    f = h5.File(filename,'r')
+    coords_dset = f["PartType0/Coordinates"]
+    coords = np.array(coords_dset)
+#translate coords by centre of box
+    header = f["Header"]
+    snap_time = header.attrs["Time"]
+    snap_time_cgs = snap_time * unit_time_cgs
+    coords[:,0] -= box_centre[0]/2.
+    coords[:,1] -= box_centre[1]/2.
+    coords[:,2] -= box_centre[2]/2.
+    radius = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2)
+    radius_cgs = radius*unit_length_cgs
+    radius_over_virial_radius = radius_cgs / r_vir_cgs
+
+#get the internal energies
+    u_dset = f["PartType0/InternalEnergy"]
+    u = np.array(u_dset)
+
+#make dimensionless
+    u /= v_c**2/(2. * (gamma - 1.))
+    r = radius_over_virial_radius
+
+    bin_edges = np.linspace(0,1,n_radial_bins + 1)
+    (hist,u_totals) = do_binning(r,u,bin_edges)
+    
+    bin_widths = 1. / n_radial_bins
+    radial_bin_mids = np.linspace(bin_widths / 2. , 1. - bin_widths / 2. , n_radial_bins) 
+    binned_u = u_totals / hist
+
+
+    plt.plot(radial_bin_mids,binned_u,'ko',label = "Numerical solution")
+    plt.plot((0,1),(1,1),label = "Analytic Solution")
+    plt.legend(loc = "lower right")
+    plt.xlabel(r"$r / r_{vir}$")
+    plt.ylabel(r"$u / (v_c^2 / (2(\gamma - 1)) $")
+    plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
+    plt.ylim((0,1))
+    plot_filename = "internal_energy_profile_%03d.png" %i
+    plt.savefig(plot_filename,format = "png")
+    plt.close()
+
+
+        
+    

examples/CoolingHalo/makeIC.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@durham.ac.uk)
+ # 
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ # 
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ # GNU General Public License for more details.
+ # 
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ # 
+ ##############################################################################
+
+import h5py
+import sys
+import numpy as np
+import math
+import random
+
+# Generates N particles in a spherically symmetric distribution with density profile ~r^(-2)
+# usage: python makeIC.py 1000: generate 1000 particles
+
+# Some constants
+
+OMEGA = 0.3 # Cosmological matter fraction at z = 0
+PARSEC_IN_CGS = 3.0856776e18
+KM_PER_SEC_IN_CGS = 1.0e5
+CONST_G_CGS = 6.672e-8
+h = 0.67777 # hubble parameter
+gamma = 5./3.
+eta = 1.2349
+
+# First set unit velocity and then the circular velocity parameter for the isothermal potential
+const_unit_velocity_in_cgs = 1.e5 #kms^-1
+
+v_c = 200.
+v_c_cgs = v_c * const_unit_velocity_in_cgs
+
+# Now we use this to get the virial mass and virial radius, which we will set to be the unit mass and radius
+
+# Find H_0, the inverse Hubble time, in cgs
+
+H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+
+# From this we can find the virial radius, the radius within which the average density of the halo is
+# 200. * the mean matter density
+
+r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
+
+# Now get the virial mass
+
+M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+
+# Now set the unit length and mass
+
+const_unit_mass_in_cgs = M_vir_cgs
+const_unit_length_in_cgs = r_vir_cgs
+
+print "UnitMass_in_cgs:     ", const_unit_mass_in_cgs 
+print "UnitLength_in_cgs:   ", const_unit_length_in_cgs
+print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
+
+#derived quantities
+const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
+print "UnitTime_in_cgs:     ", const_unit_time_in_cgs
+const_G                = ((CONST_G_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
+print 'G=', const_G
+
+# Parameters
+periodic= 1            # 1 For periodic box
+boxSize = 4.          
+G       = const_G 
+N       = int(sys.argv[1])  # Number of particles
+
+# Create the file
+filename = "CoolingHalo.hdf5"
+file = h5py.File(filename, 'w')
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs 
+grp.attrs["Unit time in cgs (U_t)"] = const_unit_length_in_cgs / const_unit_velocity_in_cgs
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+
+# Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+# set seed for random number
+np.random.seed(1234)
+
+
+# Positions
+# r^(-2) distribution corresponds to uniform distribution in radius
+radius = boxSize  * np.sqrt(3.) / 2.* np.random.rand(N) #the diagonal extent of the cube
+ctheta = -1. + 2 * np.random.rand(N)
+stheta = np.sqrt(1.-ctheta**2)
+phi    =  2 * math.pi * np.random.rand(N)
+coords      = np.zeros((N, 3))
+coords[:,0] = radius * stheta * np.cos(phi)
+coords[:,1] = radius * stheta * np.sin(phi)
+coords[:,2] = radius * ctheta
+
+#shift to centre of box
+coords += np.full((N,3),boxSize/2.)
+print "x range = (%f,%f)" %(np.min(coords[:,0]),np.max(coords[:,0]))
+print "y range = (%f,%f)" %(np.min(coords[:,1]),np.max(coords[:,1]))
+print "z range = (%f,%f)" %(np.min(coords[:,2]),np.max(coords[:,2]))
+
+print np.mean(coords[:,0])
+print np.mean(coords[:,1])
+print np.mean(coords[:,2])
+
+#now find the particles which are within the box
+
+x_coords = coords[:,0]
+y_coords = coords[:,1]
+z_coords = coords[:,2]
+
+ind = np.where(x_coords < boxSize)[0]
+x_coords = x_coords[ind]
+y_coords = y_coords[ind]
+z_coords = z_coords[ind]
+
+ind = np.where(x_coords > 0.)[0]
+x_coords = x_coords[ind]
+y_coords = y_coords[ind]
+z_coords = z_coords[ind]
+
+ind = np.where(y_coords < boxSize)[0]
+x_coords = x_coords[ind]
+y_coords = y_coords[ind]
+z_coords = z_coords[ind]
+
+ind = np.where(y_coords > 0.)[0]
+x_coords = x_coords[ind]
+y_coords = y_coords[ind]
+z_coords = z_coords[ind]
+
+ind = np.where(z_coords < boxSize)[0]
+x_coords = x_coords[ind]
+y_coords = y_coords[ind]
+z_coords = z_coords[ind]
+
+ind = np.where(z_coords > 0.)[0]
+x_coords = x_coords[ind]
+y_coords = y_coords[ind]
+z_coords = z_coords[ind]
+
+#count number of particles
+
+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
+
+radius = np.sqrt(coords_2[:,0]**2 + coords_2[:,1]**2 + coords_2[:,2]**2)
+
+#test we've done it right
+
+print "x range = (%f,%f)" %(np.min(coords_2[:,0]),np.max(coords_2[:,0]))
+print "y range = (%f,%f)" %(np.min(coords_2[:,1]),np.max(coords_2[:,1]))
+print "z range = (%f,%f)" %(np.min(coords_2[:,2]),np.max(coords_2[:,2]))
+
+print np.mean(coords_2[:,0])
+print np.mean(coords_2[:,1])
+print np.mean(coords_2[:,2])
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] =  [N ,0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [N, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["Dimension"] = 3
+
+# Particle group
+grp = file.create_group("/PartType0")
+
+ds = grp.create_dataset('Coordinates', (N, 3), 'd')
+ds[()] = coords_2
+coords_2 = np.zeros(1)
+
+# All velocities set to zero
+v = np.zeros((N,3))
+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)
+ds = grp.create_dataset('Masses', (N, ), 'f')
+ds[()] = m
+m = np.zeros(1)
+
+# Smoothing lengths
+l = (4. * np.pi * radius**2 / N)**(1./3.) #local mean inter-particle separation
+h = np.full((N, ), eta * l)
+ds = grp.create_dataset('SmoothingLength', (N,), 'f')
+ds[()] = h
+h = np.zeros(1)
+
+# Internal energies
+u = v_c**2 / (2. * (gamma - 1.)) 
+u = np.full((N, ), u)
+ds = grp.create_dataset('InternalEnergy', (N,), 'f')
+ds[()] = u
+u = np.zeros(1)
+
+# Particle IDs
+ids = 1 + np.linspace(0, N, N, endpoint=False, dtype='L')
+ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
+ds[()] = ids
+
+file.close()

examples/CoolingHalo/makeIC_random_box.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 Stefan Arridge (stefan.arridge@durham.ac.uk)
+ # 
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ # 
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ # GNU General Public License for more details.
+ # 
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ # 
+ ##############################################################################
+
+import h5py
+import sys
+import numpy as np
+import math
+import random
+
+# Generates N particles in a spherically symmetric distribution with density profile ~r^(-2)
+# usage: python makeIC.py 1000: generate 1000 particles
+
+# Some constants
+
+OMEGA = 0.3 # Cosmological matter fraction at z = 0
+PARSEC_IN_CGS = 3.0856776e18
+KM_PER_SEC_IN_CGS = 1.0e5
+CONST_G_CGS = 6.672e-8
+h = 0.67777 # hubble parameter
+gamma = 5./3.
+eta = 1.2349
+
+# First set unit velocity and then the circular velocity parameter for the isothermal potential
+const_unit_velocity_in_cgs = 1.e5 #kms^-1
+
+v_c = 200.
+v_c_cgs = v_c * const_unit_velocity_in_cgs
+
+# Now we use this to get the virial mass and virial radius, which we will set to be the unit mass and radius
+
+# Find H_0, the inverse Hubble time, in cgs
+
+H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+
+# From this we can find the virial radius, the radius within which the average density of the halo is
+# 200. * the mean matter density
+
+r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
+
+# Now get the virial mass
+
+M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+
+# Now set the unit length and mass
+
+const_unit_mass_in_cgs = M_vir_cgs
+const_unit_length_in_cgs = r_vir_cgs
+
+print "UnitMass_in_cgs:     ", const_unit_mass_in_cgs 
+print "UnitLength_in_cgs:   ", const_unit_length_in_cgs
+print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
+
+#derived quantities
+const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
+print "UnitTime_in_cgs:     ", const_unit_time_in_cgs
+const_G                = ((CONST_G_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
+print 'G=', const_G
+
+# Parameters
+periodic= 1            # 1 For periodic box
+boxSize = 4.          
+G       = const_G 
+N       = int(sys.argv[1])  # Number of particles
+
+# Create the file
+filename = "random_box.hdf5"
+file = h5py.File(filename, 'w')
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
+grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs 
+grp.attrs["Unit time in cgs (U_t)"] = const_unit_length_in_cgs / const_unit_velocity_in_cgs
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] =  [N ,0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [N, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["Dimension"] = 3
+
+# Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+# set seed for random number
+np.random.seed(1234)
+
+# Particle group
+grp = file.create_group("/PartType0")
+
+# Positions
+
+# distribute particles randomly in the box
+coords      = np.zeros((N, 3))
+coords[:,0] = boxSize * np.random.rand(N)
+coords[:,1] = boxSize * np.random.rand(N)
+coords[:,2] = boxSize * np.random.rand(N)
+#shift to centre of box
+#coords += np.full((N,3),boxSize/2.)
+print "x range = (%f,%f)" %(np.min(coords[:,0]),np.max(coords[:,0]))
+print "y range = (%f,%f)" %(np.min(coords[:,1]),np.max(coords[:,1]))
+print "z range = (%f,%f)" %(np.min(coords[:,2]),np.max(coords[:,2]))
+
+print np.mean(coords[:,0])
+print np.mean(coords[:,1])
+print np.mean(coords[:,2])
+
+ds = grp.create_dataset('Coordinates', (N, 3), 'd')
+ds[()] = coords
+coords = np.zeros(1)
+
+# All velocities set to zero
+v = np.zeros((N,3))
+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)
+ds = grp.create_dataset('Masses', (N, ), 'f')
+ds[()] = m
+m = np.zeros(1)
+
+# Smoothing lengths
+l = (boxSize**3 / N)**(1./3.) #local mean inter-particle separation
+h = np.full((N, ), eta * l)
+ds = grp.create_dataset('SmoothingLength', (N,), 'f')
+ds[()] = h
+h = np.zeros(1)
+
+# Internal energies
+u = v_c**2 / (2. * (gamma - 1.)) 
+u = np.full((N, ), u)
+ds = grp.create_dataset('InternalEnergy', (N,), 'f')
+ds[()] = u
+u = np.zeros(1)
+
+# Particle IDs
+ids = 1 + np.linspace(0, N, N, endpoint=False, dtype='L')
+ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
+ds[()] = ids
+
+file.close()

examples/CoolingHalo/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+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
+
+python radial_profile.py 2. 200 100
+
+python internal_energy_profile.py 2. 200 100
+
+python test_energy_conservation.py 2. 200 100

examples/CoolingHalo/test_energy_conservation.py

+import numpy as np
+import h5py as h5
+import matplotlib.pyplot as plt
+import sys
+
+n_snaps = 41
+
+#some constants
+OMEGA = 0.3 # Cosmological matter fraction at z = 0
+PARSEC_IN_CGS = 3.0856776e18
+KM_PER_SEC_IN_CGS = 1.0e5
+CONST_G_CGS = 6.672e-8
+h = 0.67777 # hubble parameter
+gamma = 5./3.
+eta = 1.2349
+H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+
+#read some header/parameter information from the first snapshot
+
+filename = "CoolingHalo_000.hdf5"
+f = h5.File(filename,'r')
+params = f["Parameters"]
+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_cgs = v_c * unit_velocity_cgs
+header = f["Header"]
+N = header.attrs["NumPart_Total"][0]
+box_centre = np.array(header.attrs["BoxSize"])
+
+#calculate r_vir and M_vir from v_c
+r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
+M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+
+potential_energy_array = []
+internal_energy_array = []
+kinetic_energy_array = []
+time_array_cgs = []
+
+for i in range(n_snaps):
+
+    filename = "CoolingHalo_%03d.hdf5" %i
+    f = h5.File(filename,'r')
+    coords_dset = f["PartType0/Coordinates"]
+    coords = np.array(coords_dset)
+#translate coords by centre of box
+    header = f["Header"]
+    snap_time = header.attrs["Time"]
+    snap_time_cgs = snap_time * unit_time_cgs
+    time_array_cgs = np.append(time_array_cgs,snap_time_cgs)
+    coords[:,0] -= box_centre[0]/2.
+    coords[:,1] -= box_centre[1]/2.
+    coords[:,2] -= box_centre[2]/2.
+    radius = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2)
+    radius_cgs = radius*unit_length_cgs
+    radius_over_virial_radius = radius_cgs / r_vir_cgs
+
+    r = radius_over_virial_radius
+    total_potential_energy = np.sum(v_c**2*np.log(r))
+    potential_energy_array = np.append(potential_energy_array,total_potential_energy)
+
+    vels_dset = f["PartType0/Velocities"]
+    vels = np.array(vels_dset)
+    speed_squared = vels[:,0]**2 + vels[:,1]**2 + vels[:,2]**2
+    total_kinetic_energy = 0.5 * np.sum(speed_squared)
+    kinetic_energy_array = np.append(kinetic_energy_array,total_kinetic_energy)
+
+    u_dset = f["PartType0/InternalEnergy"]
+    u = np.array(u_dset)
+    total_internal_energy = np.sum(u)
+    internal_energy_array = np.append(internal_energy_array,total_internal_energy)
+
+#put energies in units of v_c^2 and rescale by number of particles
+
+pe = potential_energy_array / (N*v_c**2)
+ke = kinetic_energy_array / (N*v_c**2)
+ie = internal_energy_array / (N*v_c**2)
+te = pe + ke + ie
+
+dyn_time_cgs = r_vir_cgs / v_c_cgs
+time_array = time_array_cgs / dyn_time_cgs
+
+plt.plot(time_array,ke,label = "Kinetic Energy")
+plt.plot(time_array,pe,label = "Potential Energy")
+plt.plot(time_array,ie,label = "Internal Energy")
+plt.plot(time_array,te,label = "Total Energy")
+plt.legend(loc = "upper right")
+plt.xlabel(r"$t / t_{dyn}$")
+plt.ylabel(r"$E / v_c^2$")
+plt.title(r"$%d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(N,v_c))
+plt.ylim((-4,2))
+#plot_filename = "density_profile_%03d.png" %i
+plt.show()
+

examples/CoolingHalo/velocity_profile.py

+import numpy as np
+import h5py as h5
+import matplotlib.pyplot as plt
+import sys
+
+def do_binning(x,y,x_bin_edges):
+
+    #x and y are arrays, where y = f(x)
+    #returns number of elements of x in each bin, and the total of the y elements corresponding to those x values
+
+    n_bins = x_bin_edges.size - 1
+    count = np.zeros(n_bins)
+    y_totals = np.zeros(n_bins)
+    
+    for i in range(n_bins):
+        ind = np.intersect1d(np.where(x > bin_edges[i])[0],np.where(x <= bin_edges[i+1])[0])
+        count[i] = ind.size
+        binned_y = y[ind]
+        y_totals[i] = np.sum(binned_y)
+
+    return(count,y_totals)
+
+
+#for the plotting
+max_r = float(sys.argv[1])
+n_radial_bins = int(sys.argv[2])
+n_snaps = int(sys.argv[3])
+
+#some constants
+OMEGA = 0.3 # Cosmological matter fraction at z = 0
+PARSEC_IN_CGS = 3.0856776e18
+KM_PER_SEC_IN_CGS = 1.0e5
+CONST_G_CGS = 6.672e-8
+CONST_m_H_CGS = 1.67e-24
+h = 0.67777 # hubble parameter
+gamma = 5./3.
+eta = 1.2349
+H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+
+#read some header/parameter information from the first snapshot
+
+filename = "CoolingHalo_000.hdf5"
+f = h5.File(filename,'r')
+params = f["Parameters"]
+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_cgs = v_c * unit_velocity_cgs
+header = f["Header"]
+N = header.attrs["NumPart_Total"][0]
+box_centre = np.array(header.attrs["BoxSize"])
+
+#calculate r_vir and M_vir from v_c
+r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
+M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+
+for i in range(n_snaps):
+
+    filename = "CoolingHalo_%03d.hdf5" %i
+    f = h5.File(filename,'r')
+    coords_dset = f["PartType0/Coordinates"]
+    coords = np.array(coords_dset)
+#translate coords by centre of box
+    header = f["Header"]
+    snap_time = header.attrs["Time"]
+    snap_time_cgs = snap_time * unit_time_cgs
+    coords[:,0] -= box_centre[0]/2.
+    coords[:,1] -= box_centre[1]/2.
+    coords[:,2] -= box_centre[2]/2.
+    radius = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2)
+    radius_cgs = radius*unit_length_cgs
+    radius_over_virial_radius = radius_cgs / r_vir_cgs
+
+#get the internal energies
+    vel_dset = f["PartType0/Velocities"]
+    vel = np.array(vel_dset)
+
+#make dimensionless
+    vel /= v_c
+    r = radius_over_virial_radius
+
+    #find radial component of velocity
+
+    v_r = np.zeros(r.size)
+    for j in range(r.size):
+        v_r[j] = -np.dot(coords[j,:],vel[j,:])/radius[j]
+
+    bin_edges = np.linspace(0,max_r,n_radial_bins + 1)
+    (hist,v_r_totals) = do_binning(r,v_r,bin_edges)
+    
+    bin_widths = bin_edges[1] - bin_edges[0]
+    radial_bin_mids = np.linspace(bin_widths / 2. , max_r - bin_widths / 2. , n_radial_bins) 
+    binned_v_r = v_r_totals / hist
+
+    #calculate cooling radius
+
+    #r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)
+
+    plt.plot(radial_bin_mids,binned_v_r,'ko',label = "Average radial velocity in shell")
+    #plt.plot((0,1),(1,1),label = "Analytic Solution")
+    #plt.plot((r_cool_over_r_vir,r_cool_over_r_vir),(0,2),'r',label = "Cooling radius")
+    plt.legend(loc = "upper right")
+    plt.xlabel(r"$r / r_{vir}$")
+    plt.ylabel(r"$v_r / v_c$")
+    plt.title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
+    plt.ylim((0,2))
+    plot_filename = "./plots/radial_velocity_profile/velocity_profile_%03d.png" %i
+    plt.savefig(plot_filename,format = "png")
+    plt.close()

examples/CoolingHaloWithSpin/README

+
+To make the initial conditions we distribute gas particles randomly in
+a cube with a side length twice that of the virial radius. The density
+profile of the gas is proportional to r^(-2) where r is the distance
+from the centre of the cube.
+
+The parameter v_rot (in makeIC.py and cooling.yml) sets the circular
+velocity of the halo, and by extension, the viral radius, viral mass,
+and the internal energy of the gas such that hydrostatic equilibrium
+is achieved.
+
+The gas is given some angular momentum about the z-axis. This is
+defined by the 'spin_lambda' variable in makeIC.py
+
+While the system is initially in hydrostatic equilibrium, the cooling
+of the gas and the non-zero angular momentum means that the halo will
+collapse into a spinning disc.
+
+To run this example, make such that the code is compiled with either
+the isothermal potential or softened isothermal potential, and
+'const_lambda' cooling, set in src/const.h. In the latter case, a
+(small) value of epsilon needs to be set in cooling.yml.  0.1 kpc
+should work well.
+
+The plotting scripts produce a plot of the density, internal energy
+and radial velocity profile for each
+snapshot. test_energy_conservation.py shows the evolution of energy
+with time. These can be used to check if the example has run properly.
+

examples/CoolingHaloWithSpin/cooling_halo.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.9885e39     # 10^6 solar masses
+  UnitLength_in_cgs:   3.0856776e21  # Kiloparsecs
+  UnitVelocity_in_cgs: 1e5           # Kilometres per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   10.    # The end time of the simulation (in internal units).
+  dt_min:     1e-7  # 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
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            CoolingHalo  # 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 for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      1.       # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+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
+  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
+
+# Cooling parameters
+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:          0.1    # Dimensionless pre-factor for the time-step condition

examples/CoolingHaloWithSpin/density_profile.py

+import numpy as np
+import h5py as h5
+import matplotlib.pyplot as plt
+import sys
+
+#for the plotting
+max_r = float(sys.argv[1]) #in units of the virial radius 
+n_radial_bins = int(sys.argv[2])
+n_snaps = int(sys.argv[3])
+
+#some constants
+OMEGA = 0.3 # Cosmological matter fraction at z = 0
+PARSEC_IN_CGS = 3.0856776e18
+KM_PER_SEC_IN_CGS = 1.0e5
+CONST_G_CGS = 6.672e-8
+CONST_m_H_CGS = 1.67e-24
+h = 0.67777 # hubble parameter
+gamma = 5./3.
+eta = 1.2349
+H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
+
+#read some header/parameter information from the first snapshot
+
+filename = "CoolingHalo_000.hdf5"
+f = h5.File(filename,'r')
+params = f["Parameters"]
+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_cgs = v_c * unit_velocity_cgs
+lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
+X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
+header = f["Header"]
+N = header.attrs["NumPart_Total"][0]
+box_centre = np.array(header.attrs["BoxSize"])
+
+#calculate r_vir and M_vir from v_c
+r_vir_cgs = v_c_cgs / (10. * H_0_cgs * np.sqrt(OMEGA))
+M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
+
+for i in range(n_snaps):
+
+    filename = "CoolingHalo_%03d.hdf5" %i
+    f = h5.File(filename,'r')
+    coords_dset = f["PartType0/Coordinates"]
+    coords = np.array(coords_dset)
+#translate coords by centre of box
+    header = f["Header"]
+    snap_time = header.attrs["Time"]
+    snap_time_cgs = snap_time * unit_time_cgs
+    coords[:,0] -= box_centre[0]/2.
+    coords[:,1] -= box_centre[1]/2.
+    coords[:,2] -= box_centre[2]/2.
+    radius = np.sqrt(coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2)
+    radius_cgs = radius*unit_length_cgs
+    radius_over_virial_radius = radius_cgs / r_vir_cgs
+
+    r = radius_over_virial_radius
+
+    bin_edges = np.linspace(0.,max_r,n_radial_bins+1)
+    bin_width = bin_edges[1] - bin_edges[0]
+    hist = np.histogram(r,bins = bin_edges)[0] # number of particles in each bin
+
+#find the mass in each radial bin
+
+    mass_dset = f["PartType0/Masses"]
+#mass of each particles should be equal
+    part_mass = np.array(mass_dset)[0]
+    part_mass_cgs = part_mass * unit_mass_cgs
+    part_mass_over_virial_mass = part_mass_cgs / M_vir_cgs 
+
+    mass_hist = hist * part_mass_over_virial_mass
+    radial_bin_mids = np.linspace(bin_width/2.,max_r - bin_width/2.,n_radial_bins)
+#volume in each radial bin
+    volume = 4.*np.pi * radial_bin_mids**2 * bin_width
+
+#now divide hist by the volume so we have a density in each bin
+
+    density = mass_hist / volume
+
+    ##read the densities
+
+    # density_dset = f["PartType0/Density"]
+    # density = np.array(density_dset)
+    # density_cgs = density * unit_mass_cgs / unit_length_cgs**3
+    # rho = density_cgs * r_vir_cgs**3 / M_vir_cgs
+
+    t = np.linspace(10./n_radial_bins,10.0,1000)
+    rho_analytic = t**(-2)/(4.*np.pi)
+
+    #calculate cooling radius
+
+    r_cool_over_r_vir = np.sqrt((2.*(gamma - 1.)*lambda_cgs*M_vir_cgs*X_H**2)/(4.*np.pi*CONST_m_H_CGS**2*v_c_cgs**2*r_vir_cgs**3))*np.sqrt(snap_time_cgs)
+
+    #initial analytic density profile
+    
+    if (i == 0):
+        r_0 = radial_bin_mids[0]
+        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.xlabel(r"$r / r_{vir}$")
+    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.xlim((radial_bin_mids[0],max_r))
+    plt.ylim((0,20))
+    plt.plot((0,max_r),(1,1))
+    #plt.xscale('log')
+    #plt.yscale('log')
+    plt.legend(loc = "upper right")
+    plot_filename = "density_profile_%03d.png" %i
+    plt.savefig(plot_filename,format = "png")
+    plt.close()
+