Commit 1608ac8d authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Updated plotting scripts for the hydrostatic halo case to be more SSH friendly...

Updated plotting scripts for the hydrostatic halo case to be more SSH friendly and to read the correct parameters in the yaml file.
parent 0fc745a2
###############################################################################
# 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 numpy as np
import h5py as h5
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use("Agg")
from pylab import *
import sys
#for the plotting
......@@ -46,7 +67,8 @@ for i in range(n_snaps):
f = h5.File(filename,'r')
coords_dset = f["PartType0/Coordinates"]
coords = np.array(coords_dset)
#translate coords by centre of box
#translate coords by centre of box
header = f["Header"]
snap_time = header.attrs["Time"]
snap_time_cgs = snap_time * unit_time_cgs
......@@ -63,58 +85,46 @@ for i in range(n_snaps):
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
#find the mass in each radial bin
mass_dset = f["PartType0/Masses"]
#mass of each particles should be equal
#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 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
#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
#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")
figure()
plot(radial_bin_mids,density/rho_analytic_init,'ko',label = "Average density of shell")
#plot(t,rho_analytic,label = "Initial analytic density profile")
xlabel(r"$r / r_{vir}$")
ylabel(r"$\rho / \rho_{init}$")
title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
xlim((radial_bin_mids[0],max_r))
ylim((0,2))
plot((0,max_r),(1,1))
legend(loc = "upper right")
plot_filename = "./plots/density_profile/density_profile_%03d.png" %i
plt.savefig(plot_filename,format = "png")
plt.close()
savefig(plot_filename,format = "png")
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
# 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 numpy as np
import h5py as h5
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use("Agg")
from pylab import *
import sys
def do_binning(x,y,x_bin_edges):
......@@ -48,8 +69,6 @@ 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
#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"])
......@@ -64,7 +83,8 @@ for i in range(n_snaps):
f = h5.File(filename,'r')
coords_dset = f["PartType0/Coordinates"]
coords = np.array(coords_dset)
#translate coords by centre of box
#translate coords by centre of box
header = f["Header"]
snap_time = header.attrs["Time"]
snap_time_cgs = snap_time * unit_time_cgs
......@@ -75,11 +95,11 @@ for i in range(n_snaps):
radius_cgs = radius*unit_length_cgs
radius_over_virial_radius = radius_cgs / r_vir_cgs
#get the internal energies
#get the internal energies
u_dset = f["PartType0/InternalEnergy"]
u = np.array(u_dset)
#make dimensionless
#make dimensionless
u /= v_c**2/(2. * (gamma - 1.))
r = radius_over_virial_radius
......@@ -90,21 +110,16 @@ for i in range(n_snaps):
radial_bin_mids = np.linspace(bin_widths / 2. , max_r - bin_widths / 2. , n_radial_bins)
binned_u = u_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_u,'ko',label = "Numerical solution")
#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 = "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,2))
figure()
plot(radial_bin_mids,binned_u,'ko',label = "Numerical solution")
legend(loc = "lower right")
xlabel(r"$r / r_{vir}$")
ylabel(r"$u / (v_c^2 / (2(\gamma - 1)) $")
title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
ylim((0,2))
plot_filename = "./plots/internal_energy/internal_energy_profile_%03d.png" %i
plt.savefig(plot_filename,format = "png")
plt.close()
savefig(plot_filename,format = "png")
close()
......
#!/bin/bash
# Generate the initial conditions if they are not present.
echo "Generating initial conditions for the isothermal potential box example..."
python makeIC.py 100000
if [ ! -e Hydrostatic.hdf5 ]
then
echo "Generating initial conditions for the isothermal potential box example..."
python makeIC.py 100000
fi
# Run for 10 dynamical times
../swift -g -s -t 2 hydrostatic.yml 2>&1 | tee output.log
../swift -g -s -t 1 hydrostatic.yml 2>&1 | tee output.log
echo "Plotting density profiles"
mkdir plots
......
###############################################################################
# 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 numpy as np
import h5py as h5
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use("Agg")
from pylab import *
import sys
n_snaps = int(sys.argv[1])
......@@ -24,7 +45,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
unit_time_cgs = unit_length_cgs / unit_velocity_cgs
v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
v_c = float(params.attrs["IsothermalPotential:vrot"])
v_c_cgs = v_c * unit_velocity_cgs
header = f["Header"]
N = header.attrs["NumPart_Total"][0]
......@@ -45,7 +66,8 @@ for i in range(n_snaps):
f = h5.File(filename,'r')
coords_dset = f["PartType0/Coordinates"]
coords = np.array(coords_dset)
#translate coords by centre of box
#translate coords by centre of box
header = f["Header"]
snap_time = header.attrs["Time"]
snap_time_cgs = snap_time * unit_time_cgs
......@@ -73,7 +95,6 @@ for i in range(n_snaps):
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)
......@@ -82,14 +103,15 @@ 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 = "lower 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((-2,2))
plt.savefig("energy_conservation.png",format = 'png')
figure()
plot(time_array,ke,label = "Kinetic Energy")
plot(time_array,pe,label = "Potential Energy")
plot(time_array,ie,label = "Internal Energy")
plot(time_array,te,label = "Total Energy")
legend(loc = "lower right")
xlabel(r"$t / t_{dyn}$")
ylabel(r"$E / v_c^2$")
title(r"$%d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(N,v_c))
ylim((-2,2))
savefig("energy_conservation.png",format = 'png')
###############################################################################
# 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 numpy as np
import h5py as h5
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use("Agg")
from pylab import *
import sys
def do_binning(x,y,x_bin_edges):
......@@ -62,7 +83,8 @@ for i in range(n_snaps):
f = h5.File(filename,'r')
coords_dset = f["PartType0/Coordinates"]
coords = np.array(coords_dset)
#translate coords by centre of box
#translate coords by centre of box
header = f["Header"]
snap_time = header.attrs["Time"]
snap_time_cgs = snap_time * unit_time_cgs
......@@ -73,16 +95,15 @@ for i in range(n_snaps):
radius_cgs = radius*unit_length_cgs
radius_over_virial_radius = radius_cgs / r_vir_cgs
#get the internal energies
#get the internal energies
vel_dset = f["PartType0/Velocities"]
vel = np.array(vel_dset)
#make dimensionless
#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]
......@@ -94,18 +115,13 @@ for i in range(n_snaps):
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))
figure()
plot(radial_bin_mids,binned_v_r,'ko',label = "Average radial velocity in shell")
legend(loc = "upper right")
xlabel(r"$r / r_{vir}$")
ylabel(r"$v_r / v_c$")
title(r"$\mathrm{Time}= %.3g \, s \, , \, %d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %(snap_time_cgs,N,v_c))
ylim((-1,1))
plot_filename = "./plots/radial_velocity_profile/velocity_profile_%03d.png" %i
plt.savefig(plot_filename,format = "png")
plt.close()
savefig(plot_filename,format = "png")
close()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment