Commit 3334941d authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' into no_specific_counters

parents 04416705 183c658f
# 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:
......@@ -36,9 +36,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
......@@ -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()
......
......@@ -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.
......
......@@ -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_fixdt -s -C -t 16 coolingBox.yml
#-C 2>&1 | tee output.log
python energy_plot.py 0
#python test_energy_conservation.py 0
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()
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.
# 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
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()
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()
###############################################################################
# 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