Added example of a halo initially in Hydrostatic equilibrium which collapses due to cooling

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.
+  max_smoothing_length:  1.      # Maximal smoothing length allowed (in internal units).
+
+# 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/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/radial_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/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 10
+
+# python internal_energy_profile.py 10
+
+# python test_energy_conservation.py

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()
+