Commit 46865d44 authored by Stefan Arridge's avatar Stefan Arridge
Browse files

Fixed hydrostatic halo example. Should work properly now

parent 870546fc
......@@ -9,7 +9,7 @@ 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).
time_end: 10.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).
......
......@@ -21,7 +21,7 @@ def do_binning(x,y,x_bin_edges):
return(count,y_totals)
n_snaps = 11
n_snaps = 100
#for the plotting
n_radial_bins = int(sys.argv[1])
......
......@@ -90,17 +90,6 @@ grp.attrs["Unit time in cgs (U_t)"] = const_unit_length_in_cgs / const_unit_velo
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")
......@@ -109,12 +98,10 @@ grp.attrs["PeriodicBoundariesOn"] = periodic
# set seed for random number
np.random.seed(1234)
# Particle group
grp = file.create_group("/PartType0")
# Positions
# r^(-2) distribution corresponds to uniform distribution in radius
radius = np.random.rand(N)
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)
......@@ -133,9 +120,84 @@ 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
coords = np.zeros(1)
ds[()] = coords_2
coords_2 = np.zeros(1)
# All velocities set to zero
v = np.zeros((N,3))
......
......@@ -6,7 +6,7 @@ import sys
n_snaps = 11
#for the plotting
n_radial_bins = int(sys.argv[1])
#n_radial_bins = int(sys.argv[1])
#some constants
OMEGA = 0.3 # Cosmological matter fraction at z = 0
......@@ -56,36 +56,45 @@ for i in range(n_snaps):
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
# 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
# #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_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
# 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
# #now divide hist by the volume so we have a density in each bin
density = mass_hist / volume
# density = mass_hist / volume
t = np.linspace(0.01,1.0,1000)
# 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(radial_bin_mids,density,'ko',label = "Numerical solution")
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.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()
......
......@@ -2,7 +2,7 @@
# Generate the initial conditions if they are not present.
echo "Generating initial conditions for the isothermal potential box example..."
python makeIC.py 10000
python makeIC_full_box.py 10000
../swift -g -s -t 16 hydrostatic.yml 2>&1 | tee output.log
......
......@@ -3,7 +3,7 @@ import h5py as h5
import matplotlib.pyplot as plt
import sys
n_snaps = 101
n_snaps = 1000
#some constants
OMEGA = 0.3 # Cosmological matter fraction at z = 0
......@@ -63,13 +63,13 @@ for i in range(n_snaps):
vels_dset = f["PartType0/Velocities"]
vels = np.array(vels_dset)
speed_squared = coords[:,0]**2 + coords[:,1]**2 + coords[:,2]**2
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 = 0.5 * np.sum(u)
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
......@@ -79,11 +79,6 @@ ke = kinetic_energy_array / (N*v_c**2)
ie = internal_energy_array / (N*v_c**2)
te = pe + ke + ie
print pe
print ke
print ie
print te
dyn_time_cgs = r_vir_cgs / v_c_cgs
time_array = time_array_cgs / dyn_time_cgs
......@@ -98,3 +93,4 @@ plt.title(r"$%d \, \, \mathrm{particles} \,,\, v_c = %.1f \, \mathrm{km / s}$" %
plt.ylim((-2,2))
#plot_filename = "density_profile_%03d.png" %i
plt.show()
......@@ -94,7 +94,8 @@
//#define EXTERNAL_POTENTIAL_NONE
//#define EXTERNAL_POTENTIAL_POINTMASS
#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
#define EXTERNAL_POTENTIAL_SOFTENED_ISOTHERMAL_POTENTIAL
//#define EXTERNAL_POTENTIAL_DISC_PATCH
/* Source terms */
......
......@@ -37,6 +37,8 @@
#include "./potential/point_mass/potential.h"
#elif defined(EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL)
#include "./potential/isothermal/potential.h"
#elif defined(EXTERNAL_POTENTIAL_SOFTENED_ISOTHERMAL_POTENTIAL)
#include "./potential/softened_isothermal/potential.h"
#elif defined(EXTERNAL_POTENTIAL_DISC_PATCH)
#include "./potential/disc_patch/potential.h"
#else
......
......@@ -112,7 +112,7 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
g->a_grav[0] = term * dx;
g->a_grav[1] = term * dy;
g->a_grav[2] = term * dz;
g->a_grav[2] = term * dz;
}
/**
......
......@@ -39,7 +39,6 @@
/* This object's header. */
#include "task.h"
/* Local headers. */
#include "atomic.h"
#include "error.h"
......
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