Skip to content
Snippets Groups Projects
Commit eb044806 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke Committed by Matthieu Schaller
Browse files

Some more changes to Zeldovich test parameters and scripts.

parent e51b15cd
Branches
Tags
2 merge requests!566Periodic gravity calculation,!565Mesh force task
...@@ -26,8 +26,8 @@ from numpy import * ...@@ -26,8 +26,8 @@ from numpy import *
gamma = 5./3. # Gas adiabatic index gamma = 5./3. # Gas adiabatic index
numPart_1D = 32 # Number of particles numPart_1D = 32 # Number of particles
lambda_i = 1.975e24 # h^-1 m (= 64 h^-1 Mpc) lambda_i = 1.975e24 # h^-1 m (= 64 h^-1 Mpc)
x_min = 0. x_min = -0.5 * lambda_i
x_max = lambda_i x_max = 0.5 * lambda_i
rho_0 = 1.8788e-26 # h^2 kg m^-3 rho_0 = 1.8788e-26 # h^2 kg m^-3
T_i = 100. # K T_i = 100. # K
H_0 = 3.24077929e-18 # s^-1 H_0 = 3.24077929e-18 # s^-1
...@@ -62,6 +62,7 @@ boxSize = x_max - x_min ...@@ -62,6 +62,7 @@ boxSize = x_max - x_min
delta_x = boxSize / numPart_1D delta_x = boxSize / numPart_1D
# Get the particle mass # Get the particle mass
a_i = 1. / (1. + z_i)
m_i = boxSize**3 * rho_0 / numPart m_i = boxSize**3 * rho_0 / numPart
# Build the arrays # Build the arrays
...@@ -78,9 +79,9 @@ for i in range(numPart_1D): ...@@ -78,9 +79,9 @@ for i in range(numPart_1D):
for k in range(numPart_1D): for k in range(numPart_1D):
index = i * numPart_1D**2 + j * numPart_1D + k index = i * numPart_1D**2 + j * numPart_1D + k
q = x_min + (i + 0.5) * delta_x q = x_min + (i + 0.5) * delta_x
coords[index,0] = q - zfac * sin(k_i * q) / k_i coords[index,0] = q - zfac * sin(k_i * q) / k_i - x_min
coords[index,1] = x_min + (j + 0.5) * delta_x coords[index,1] = (j + 0.5) * delta_x
coords[index,2] = x_min + (k + 0.5) * delta_x coords[index,2] = (k + 0.5) * delta_x
T = T_i * (1. / (1. - zfac * cos(k_i * q)))**(2. / 3.) T = T_i * (1. / (1. - zfac * cos(k_i * q)))**(2. / 3.)
u[index] = k_in_J_K * T / (gamma - 1.) / mH_in_kg u[index] = k_in_J_K * T / (gamma - 1.) / mH_in_kg
h[index] = 1.2348 * delta_x h[index] = 1.2348 * delta_x
...@@ -135,3 +136,8 @@ grp.create_dataset('InternalEnergy', data=u, dtype='f') ...@@ -135,3 +136,8 @@ grp.create_dataset('InternalEnergy', data=u, dtype='f')
grp.create_dataset('ParticleIDs', data=ids, dtype='L') grp.create_dataset('ParticleIDs', data=ids, dtype='L')
file.close() file.close()
import pylab as pl
pl.plot(coords[:,0], v[:,0], "k.")
pl.show()
...@@ -21,14 +21,11 @@ ...@@ -21,14 +21,11 @@
# the simulation result # the simulation result
# Parameters # Parameters
gas_gamma = 5./3. # Polytropic index
lambda_i = 1.975e24 # Wavelength of the density perturbation (in h^-1 m)
rho_0 = 1.8788e-26 # Critical density of the Universe (in h^2 kg m^-3)
T_i = 100. # Initial temperature of the gas (in K) T_i = 100. # Initial temperature of the gas (in K)
H_0 = 3.24077929e-18 # Assumed value of the Hubble constant (in s^-1)
z_c = 1. # Redshift of caustic formation (non-linear collapse) z_c = 1. # Redshift of caustic formation (non-linear collapse)
z_i = 100. # Initial redshift z_i = 100. # Initial redshift
# Physical constants needed for internal energy to temperature conversion
k_in_J_K = 1.38064852e-23 k_in_J_K = 1.38064852e-23
mH_in_kg = 1.6737236e-27 mH_in_kg = 1.6737236e-27
...@@ -45,7 +42,7 @@ params = {'axes.labelsize': 10, ...@@ -45,7 +42,7 @@ params = {'axes.labelsize': 10,
'xtick.labelsize': 10, 'xtick.labelsize': 10,
'ytick.labelsize': 10, 'ytick.labelsize': 10,
'text.usetex': True, 'text.usetex': True,
'figure.figsize' : (9.90,6.45), 'figure.figsize' : (6.45,6.45),
'figure.subplot.left' : 0.045, 'figure.subplot.left' : 0.045,
'figure.subplot.right' : 0.99, 'figure.subplot.right' : 0.99,
'figure.subplot.bottom' : 0.05, 'figure.subplot.bottom' : 0.05,
...@@ -73,14 +70,26 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"] ...@@ -73,14 +70,26 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
eta = sim["/HydroScheme"].attrs["Kernel eta"] eta = sim["/HydroScheme"].attrs["Kernel eta"]
git = sim["Code"].attrs["Git Revision"] git = sim["Code"].attrs["Git Revision"]
# Cosmological parameters
H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
x = sim["/PartType0/Coordinates"][:,0] x = sim["/PartType0/Coordinates"][:,0]
v = sim["/PartType0/Velocities"][:,0] v = sim["/PartType0/Velocities"][:,0]
u = sim["/PartType0/InternalEnergy"][:] u = sim["/PartType0/InternalEnergy"][:]
S = sim["/PartType0/Entropy"][:] S = sim["/PartType0/Entropy"][:]
P = sim["/PartType0/Pressure"][:] P = sim["/PartType0/Pressure"][:]
rho = sim["/PartType0/Density"][:] rho = sim["/PartType0/Density"][:]
m = sim["/PartType0/Masses"][:]
phi = sim["/PartType0/Potential"][:]
x -= 0.5 * boxSize
x_s = linspace(0., lambda_i, 1000) # Derived parameters
rho_0 = m.sum() / (boxSize**3) # critical density of the box
lambda_i = boxSize # wavelength of the perturbation
x_s = linspace(-0.5 * lambda_i, 0.5 * lambda_i, 1000)
k_i = 2. * pi / lambda_i k_i = 2. * pi / lambda_i
zfac = (1. + z_c) / (1. + redshift) zfac = (1. + z_c) / (1. + redshift)
rho_s = rho_0 / (1 - zfac * cos(k_i * x_s)) rho_s = rho_0 / (1 - zfac * cos(k_i * x_s))
...@@ -99,56 +108,38 @@ unit_length_in_si = 0.01 * unit_length_in_cgs ...@@ -99,56 +108,38 @@ unit_length_in_si = 0.01 * unit_length_in_cgs
unit_mass_in_si = 0.001 * unit_mass_in_cgs unit_mass_in_si = 0.001 * unit_mass_in_cgs
unit_time_in_si = unit_time_in_cgs unit_time_in_si = unit_time_in_cgs
x_s /= unit_length_in_si
rho_s *= (unit_length_in_si**3 / unit_mass_in_si)
v_s *= (unit_time_in_si / unit_length_in_si)
# Plot the interesting quantities # Plot the interesting quantities
figure() figure()
# Velocity profile -------------------------------- # Velocity profile --------------------------------
subplot(231) subplot(221)
plot(x, v, '.', color='r', ms=4.0) plot(x, v, '.', color='r', ms=4.0)
plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2) plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", labelpad=0) xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Velocity}}~v_x$", labelpad=0) ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
# Density profile -------------------------------- # Density profile --------------------------------
subplot(232) subplot(222)
plot(x, rho, '.', color='r', ms=4.0) plot(x, phi, '.', color='r', ms=4.0)
plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2) #plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", labelpad=0) xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Density}}~\\rho$", labelpad=0) ylabel("${\\rm{Density}}~\\rho$", labelpad=0)
# Pressure profile --------------------------------
subplot(233)
plot(x, P, '.', color='r', ms=4.0)
plot(x_s, P_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Pressure}}~P$", labelpad=0)
# Internal energy profile ------------------------- # Internal energy profile -------------------------
subplot(234) subplot(223)
#plot(x, u, '.', color='r', ms=4.0) #plot(x, u, '.', color='r', ms=4.0)
#plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2) #plot(x_s, u_s, '--', color='k', alpha=0.8, lw=1.2)
u *= (unit_length_in_si**2 / unit_time_in_si**2) u *= (unit_length_in_si**2 / unit_time_in_si**2)
u /= a**(3 * (gas_gamma - 1.)) u /= a**(3 * (gas_gamma - 1.))
T = (gas_gamma - 1.) * u * mH_in_kg / k_in_J_K T = (gas_gamma - 1.) * u * mH_in_kg / k_in_J_K
print T.mean() print "z = {0:.2f}, T_avg = {1:.2f}".format(redshift, T.mean())
plot(x, T, '.', color='r', ms=4.0) plot(x, T, '.', color='r', ms=4.0)
plot(x_s, T_s, '--', color='k', alpha=0.8, lw=1.2) plot(x_s, T_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", labelpad=0) xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0) ylabel("${\\rm{Internal~Energy}}~u$", labelpad=0)
# Entropy profile ---------------------------------
subplot(235)
plot(x, S, '.', color='r', ms=4.0)
plot(x_s, s_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Entropy}}~S$", labelpad=0)
# Information ------------------------------------- # Information -------------------------------------
subplot(236, frameon=False) subplot(224, frameon=False)
text(-0.49, 0.9, "Zeldovich pancake with $\\gamma=%.3f$ in 1D at $t=%.2f$"%(gas_gamma,time), fontsize=10) text(-0.49, 0.9, "Zeldovich pancake with $\\gamma=%.3f$ in 1D at $t=%.2f$"%(gas_gamma,time), fontsize=10)
text(-0.49, 0.8, "$z={0:.2f}$".format(redshift)) text(-0.49, 0.8, "$z={0:.2f}$".format(redshift))
......
...@@ -17,8 +17,8 @@ TimeIntegration: ...@@ -17,8 +17,8 @@ TimeIntegration:
Snapshots: Snapshots:
basename: zeldovichPancake # Common part of the name of output files basename: zeldovichPancake # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units) time_first: 0. # Time of the first output (in internal units)
delta_time: 10 # Time difference between consecutive outputs (in internal units) delta_time: 1.04 # Time difference between consecutive outputs (in internal units)
scale_factor_first: 0.01 scale_factor_first: 0.00991
# Parameters governing the conserved quantities statistics # Parameters governing the conserved quantities statistics
Statistics: Statistics:
...@@ -36,7 +36,7 @@ InitialConditions: ...@@ -36,7 +36,7 @@ InitialConditions:
Cosmology: Cosmology:
Omega_m: 1. Omega_m: 1.
Omega_lambda: 0. Omega_lambda: 0.
Omega_b: 0. Omega_b: 1.
h: 1. h: 1.
a_begin: 0.00990099 a_begin: 0.00990099
a_end: 1.0 a_end: 1.0
...@@ -44,5 +44,5 @@ Cosmology: ...@@ -44,5 +44,5 @@ Cosmology:
Gravity: Gravity:
eta: 0.025 eta: 0.025
theta: 0.9 theta: 0.9
comoving_softening: 0.001 comoving_softening: 0.0001
max_physical_softening: 0.001 max_physical_softening: 0.0001
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment