Skip to content
Snippets Groups Projects
Commit cfd6f2c0 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'cosmo_glass' into 'master'

Use a gravity glass file to generate the ICs of the constant cosmological volume

See merge request !602
parents 1311f784 6c5d2314
No related branches found
No related tags found
1 merge request!602Use a gravity glass file to generate the ICs of the constant cosmological volume
This test is a small cosmological volume with constant density and internal energy.
The ICs are generated from a glass file to minimize the build-up of peculiar velocities
over time.
The solution script plots the expected solution both in comoving and physical frames.
...@@ -17,7 +17,7 @@ Cosmology: ...@@ -17,7 +17,7 @@ Cosmology:
# Parameters governing the time integration # Parameters governing the time integration
TimeIntegration: TimeIntegration:
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units). dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
dt_max: 5e-3 # The maximal time-step size of the simulation (in internal units). dt_max: 2e-3 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots # Parameters governing the snapshots
Snapshots: Snapshots:
......
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/gravity_glassCube_32.hdf5
...@@ -25,7 +25,7 @@ T_i = 100. # Initial temperature of the gas (in K) ...@@ -25,7 +25,7 @@ T_i = 100. # Initial temperature of the gas (in K)
z_i = 100. # Initial redshift z_i = 100. # Initial redshift
gamma = 5./3. # Gas adiabatic index gamma = 5./3. # Gas adiabatic index
numPart_1D = 32 numPart_1D = 32
#glassFile = "glassCube_32.hdf5" glassFile = "gravity_glassCube_32.hdf5"
fileName = "constantBox.hdf5" fileName = "constantBox.hdf5"
...@@ -56,17 +56,16 @@ unit_u_in_si = unit_v_in_si**2 ...@@ -56,17 +56,16 @@ unit_u_in_si = unit_v_in_si**2
#--------------------------------------------------- #---------------------------------------------------
# Read the glass file # Read the glass file
#glass = h5py.File(glassFile, "r" ) glass = h5py.File(glassFile, "r" )
# Read particle positions and h from the glass # Read particle positions and h from the glass
#pos = glass["/PartType0/Coordinates"][:,:] pos = glass["/PartType1/Coordinates"][:,:]
#h = glass["/PartType0/SmoothingLength"][:] * 0.3 glass.close()
#glass.close()
# Total number of particles # Total number of particles
#numPart = size(h) numPart = size(pos)/3
#if numPart != numPart_1D**3: if numPart != numPart_1D**3:
# print "Non-matching glass file" print "Non-matching glass file"
numPart = numPart_1D**3 numPart = numPart_1D**3
# Set box size and interparticle distance # Set box size and interparticle distance
...@@ -78,9 +77,7 @@ a_i = 1. / (1. + z_i) ...@@ -78,9 +77,7 @@ 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
#pos *= boxSize pos *= boxSize
#h *= boxSize
coords = zeros((numPart, 3))
v = zeros((numPart, 3)) v = zeros((numPart, 3))
ids = linspace(1, numPart, numPart) ids = linspace(1, numPart, numPart)
m = zeros(numPart) m = zeros(numPart)
...@@ -92,9 +89,9 @@ for i in range(numPart_1D): ...@@ -92,9 +89,9 @@ for i in range(numPart_1D):
for j in range(numPart_1D): for j 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
coords[index,0] = (i + 0.5) * delta_x #coords[index,0] = (i + 0.5) * delta_x
coords[index,1] = (j + 0.5) * delta_x #coords[index,1] = (j + 0.5) * delta_x
coords[index,2] = (k + 0.5) * delta_x #coords[index,2] = (k + 0.5) * delta_x
u[index] = kB_in_SI * T_i / (gamma - 1.) / mH_in_kg u[index] = kB_in_SI * T_i / (gamma - 1.) / mH_in_kg
h[index] = 1.2348 * delta_x h[index] = 1.2348 * delta_x
m[index] = m_i m[index] = m_i
...@@ -103,7 +100,7 @@ for i in range(numPart_1D): ...@@ -103,7 +100,7 @@ for i in range(numPart_1D):
v[index,2] = 0. v[index,2] = 0.
# Unit conversion # Unit conversion
coords /= unit_l_in_si pos /= unit_l_in_si
v /= unit_v_in_si v /= unit_v_in_si
m /= unit_m_in_si m /= unit_m_in_si
h /= unit_l_in_si h /= unit_l_in_si
...@@ -140,7 +137,7 @@ grp.attrs["Unit temperature in cgs (U_T)"] = 1. ...@@ -140,7 +137,7 @@ grp.attrs["Unit temperature in cgs (U_T)"] = 1.
#Particle group #Particle group
grp = file.create_group("/PartType0") grp = file.create_group("/PartType0")
grp.create_dataset('Coordinates', data=coords, dtype='d', compression="gzip", shuffle=True) grp.create_dataset('Coordinates', data=pos, dtype='d', compression="gzip", shuffle=True)
grp.create_dataset('Velocities', data=v, dtype='f',compression="gzip", shuffle=True) grp.create_dataset('Velocities', data=v, dtype='f',compression="gzip", shuffle=True)
grp.create_dataset('Masses', data=m, dtype='f', compression="gzip", shuffle=True) grp.create_dataset('Masses', data=m, dtype='f', compression="gzip", shuffle=True)
grp.create_dataset('SmoothingLength', data=h, dtype='f', compression="gzip", shuffle=True) grp.create_dataset('SmoothingLength', data=h, dtype='f', compression="gzip", shuffle=True)
......
#!/bin/bash #!/bin/bash
# Generate the initial conditions if they are not present. # Generate the initial conditions if they are not present.
if [ ! -e gravity_glassCube_32.hdf5 ]
then
echo "Fetching initial grvity glass file for the constant cosmological box example..."
./getGlass.sh
fi
if [ ! -e constantBox.hdf5 ] if [ ! -e constantBox.hdf5 ]
then then
echo "Generating initial conditions for the uniform cosmo box example..." echo "Generating initial conditions for the uniform cosmo box example..."
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment