From eb0448064b26048bc4cce94049a231da3dbbf762 Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@gmail.com>
Date: Wed, 23 May 2018 13:48:24 +0100
Subject: [PATCH] Some more changes to Zeldovich test parameters and scripts.
---
examples/ZeldovichPancake_3D/makeIC.py | 16 ++++--
examples/ZeldovichPancake_3D/plotSolution.py | 53 ++++++++-----------
.../ZeldovichPancake_3D/zeldovichPancake.yml | 10 ++--
3 files changed, 38 insertions(+), 41 deletions(-)
diff --git a/examples/ZeldovichPancake_3D/makeIC.py b/examples/ZeldovichPancake_3D/makeIC.py
index a191e3a0d7..4e660b02d3 100644
--- a/examples/ZeldovichPancake_3D/makeIC.py
+++ b/examples/ZeldovichPancake_3D/makeIC.py
@@ -26,8 +26,8 @@ from numpy import *
gamma = 5./3. # Gas adiabatic index
numPart_1D = 32 # Number of particles
lambda_i = 1.975e24 # h^-1 m (= 64 h^-1 Mpc)
-x_min = 0.
-x_max = lambda_i
+x_min = -0.5 * lambda_i
+x_max = 0.5 * lambda_i
rho_0 = 1.8788e-26 # h^2 kg m^-3
T_i = 100. # K
H_0 = 3.24077929e-18 # s^-1
@@ -62,6 +62,7 @@ boxSize = x_max - x_min
delta_x = boxSize / numPart_1D
# Get the particle mass
+a_i = 1. / (1. + z_i)
m_i = boxSize**3 * rho_0 / numPart
# Build the arrays
@@ -78,9 +79,9 @@ for i in range(numPart_1D):
for k in range(numPart_1D):
index = i * numPart_1D**2 + j * numPart_1D + k
q = x_min + (i + 0.5) * delta_x
- coords[index,0] = q - zfac * sin(k_i * q) / k_i
- coords[index,1] = x_min + (j + 0.5) * delta_x
- coords[index,2] = x_min + (k + 0.5) * delta_x
+ coords[index,0] = q - zfac * sin(k_i * q) / k_i - x_min
+ coords[index,1] = (j + 0.5) * delta_x
+ coords[index,2] = (k + 0.5) * delta_x
T = T_i * (1. / (1. - zfac * cos(k_i * q)))**(2. / 3.)
u[index] = k_in_J_K * T / (gamma - 1.) / mH_in_kg
h[index] = 1.2348 * delta_x
@@ -135,3 +136,8 @@ grp.create_dataset('InternalEnergy', data=u, dtype='f')
grp.create_dataset('ParticleIDs', data=ids, dtype='L')
file.close()
+
+import pylab as pl
+
+pl.plot(coords[:,0], v[:,0], "k.")
+pl.show()
diff --git a/examples/ZeldovichPancake_3D/plotSolution.py b/examples/ZeldovichPancake_3D/plotSolution.py
index 0335483213..fff1fa06fb 100644
--- a/examples/ZeldovichPancake_3D/plotSolution.py
+++ b/examples/ZeldovichPancake_3D/plotSolution.py
@@ -21,14 +21,11 @@
# the simulation result
# 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)
-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_i = 100. # Initial redshift
+# Physical constants needed for internal energy to temperature conversion
k_in_J_K = 1.38064852e-23
mH_in_kg = 1.6737236e-27
@@ -45,7 +42,7 @@ params = {'axes.labelsize': 10,
'xtick.labelsize': 10,
'ytick.labelsize': 10,
'text.usetex': True,
- 'figure.figsize' : (9.90,6.45),
+ 'figure.figsize' : (6.45,6.45),
'figure.subplot.left' : 0.045,
'figure.subplot.right' : 0.99,
'figure.subplot.bottom' : 0.05,
@@ -73,14 +70,26 @@ neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"]
eta = sim["/HydroScheme"].attrs["Kernel eta"]
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]
v = sim["/PartType0/Velocities"][:,0]
u = sim["/PartType0/InternalEnergy"][:]
S = sim["/PartType0/Entropy"][:]
P = sim["/PartType0/Pressure"][:]
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
zfac = (1. + z_c) / (1. + redshift)
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
unit_mass_in_si = 0.001 * unit_mass_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
figure()
# Velocity profile --------------------------------
-subplot(231)
+subplot(221)
plot(x, v, '.', color='r', ms=4.0)
plot(x_s, v_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", labelpad=0)
ylabel("${\\rm{Velocity}}~v_x$", labelpad=0)
# Density profile --------------------------------
-subplot(232)
-plot(x, rho, '.', color='r', ms=4.0)
-plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
+subplot(222)
+plot(x, phi, '.', color='r', ms=4.0)
+#plot(x_s, rho_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", 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 -------------------------
-subplot(234)
+subplot(223)
#plot(x, u, '.', color='r', ms=4.0)
#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 /= a**(3 * (gas_gamma - 1.))
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_s, T_s, '--', color='k', alpha=0.8, lw=1.2)
xlabel("${\\rm{Position}}~x$", 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 -------------------------------------
-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.8, "$z={0:.2f}$".format(redshift))
diff --git a/examples/ZeldovichPancake_3D/zeldovichPancake.yml b/examples/ZeldovichPancake_3D/zeldovichPancake.yml
index dd70eb7154..e386a7fc09 100644
--- a/examples/ZeldovichPancake_3D/zeldovichPancake.yml
+++ b/examples/ZeldovichPancake_3D/zeldovichPancake.yml
@@ -17,8 +17,8 @@ TimeIntegration:
Snapshots:
basename: zeldovichPancake # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
- delta_time: 10 # Time difference between consecutive outputs (in internal units)
- scale_factor_first: 0.01
+ delta_time: 1.04 # Time difference between consecutive outputs (in internal units)
+ scale_factor_first: 0.00991
# Parameters governing the conserved quantities statistics
Statistics:
@@ -36,7 +36,7 @@ InitialConditions:
Cosmology:
Omega_m: 1.
Omega_lambda: 0.
- Omega_b: 0.
+ Omega_b: 1.
h: 1.
a_begin: 0.00990099
a_end: 1.0
@@ -44,5 +44,5 @@ Cosmology:
Gravity:
eta: 0.025
theta: 0.9
- comoving_softening: 0.001
- max_physical_softening: 0.001
+ comoving_softening: 0.0001
+ max_physical_softening: 0.0001
--
GitLab