Commit 6e2c2eac authored by Stefan Arridge's avatar Stefan Arridge
Browse files

Added 'C' cooling flag as an option

parent fdf6cc98
......@@ -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: 8. # The end time of the simulation (in internal units).
time_end: 10. # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-1 # The maximal time-step size of the simulation (in internal units).
......@@ -21,7 +21,7 @@ Statistics:
Snapshots:
basename: Isothermal # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.02 # Time difference between consecutive outputs (in internal units)
delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
# Parameters for the hydrodynamics scheme
SPH:
......@@ -33,15 +33,15 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: Isothermal.hdf5 # The file to read
shift_x: 100. # A shift to apply to all particles read from the ICs (in internal units).
shift_y: 100.
shift_z: 100.
shift_x: 200. # A shift to apply to all particles read from the ICs (in internal units).
shift_y: 200.
shift_z: 200.
# External potential parameters
IsothermalPotential:
position_x: 100. # location of centre of isothermal potential in internal units
position_y: 100.
position_z: 100.
position_x: 200. # location of centre of isothermal potential in internal units
position_y: 200.
position_z: 200.
vrot: 200. # rotation speed of isothermal potential in internal units
timestep_mult: 0.03 # controls time step
......@@ -25,8 +25,8 @@ import math
import random
# Generates N particles in a spherical distribution centred on [0,0,0], to be moved in an isothermal potential
# usage: python makeIC.py 1000 0 : generate 1000 particles on circular orbits
# python makeIC.py 1000 1 : generate 1000 particles with Lz/L uniform in [0,1]
# usage: python makeIC.py 1000 0 1: generate 1000 particles on circular orbits
# python makeIC.py 1000 1 1: generate 1000 particles with Lz/L uniform in [0,1]
# all particles move in the xy plane, and start at y=0
# physical constants in cgs
......@@ -64,6 +64,7 @@ G = const_G
N = int(sys.argv[1]) # Number of particles
icirc = int(sys.argv[2]) # if = 0, all particles are on circular orbits, if = 1, Lz/Lcirc uniform in ]0,1[
three_d = int(sys.argv[3]) # if = 0, particles distributed in 3 dimensions, if 1, all along one straight line
L = N**(1./3.)
# these are not used but necessary for I/O
......@@ -119,23 +120,44 @@ ctheta = -1. + 2 * numpy.random.rand(N)
stheta = numpy.sqrt(1.-ctheta**2)
phi = 2 * math.pi * numpy.random.rand(N)
r = numpy.zeros((numPart, 3))
#r[:,0] = radius * stheta * numpy.cos(phi)
#r[:,1] = radius * stheta * numpy.sin(phi)
#r[:,2] = radius * ctheta
r[:,0] = radius
if (three_d == 0):
r[:,0] = radius * stheta * numpy.cos(phi)
r[:,1] = radius * stheta * numpy.sin(phi)
r[:,2] = radius * ctheta
else:
r[:,0] = radius
#
speed = vrot
v = numpy.zeros((numPart, 3))
omega = speed / radius
period = 2.*math.pi/omega
print 'period = minimum = ',min(period), ' maximum = ',max(period)
omegav = omega
if (icirc != 0):
omegav = omega * numpy.random.rand(N)
v[:,0] = -omegav * r[:,1]
v[:,1] = omegav * r[:,0]
if(three_d == 0):
#make N random angular momentum vectors
omega = numpy.random.rand(N,3)
#cross with r to get velocities perpendicular to displacement vector
v = numpy.cross(omega,r)
#rescale v
mag_v = numpy.zeros(N)
for i in range(N):
mag_v[i] = numpy.sqrt(v[i,0]**2 + v[i,1]**2 + v[i,2]**2)
v[i,:] /= mag_v[i]
v *= speed
if (icirc != 0):
for i in range(N):
v[i,:] *= numpy.random.rand()
else:
omega = speed / radius
period = 2.*math.pi/omega
print 'period = minimum = ',min(period), ' maximum = ',max(period)
omegav = omega
if (icirc != 0):
omegav = omega * numpy.random.rand(N)
v[:,0] = -omegav * r[:,1]
v[:,1] = omegav * r[:,0]
ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
......
import numpy as np
import matplotlib.pyplot as plt
import h5py as h5
n_snaps = 102
directory = '../'
snap_name_base = "Isothermal_"
plot_dir = "random_orbits_1d"
plot_filename_base = "snap"
v_rot = 200.
r_0 = 100.
e_kin_array = []
e_pot_array = []
e_total_array = []
snap_numbers = np.linspace(0,n_snaps)
for i in range(n_snaps):
snap_number = "%03d" %i
filename = directory + snap_name_base + snap_number + ".hdf5"
f = h5.File(filename)
header = f["Header"]
n_parts = header.attrs["NumPart_ThisFile"][1]
pos = np.array(f["PartType1/Coordinates"])
x = pos[:,0]
y = pos[:,1]
z = pos[:,2]
r = np.sqrt(x**2 + y**2 + z**2)
e_pot_particles = v_rot**2 * np.log(r/r_0)
e_pot_total = np.sum(e_pot_particles)
e_pot_array = np.append(e_pot_array,e_pot_total)
v = np.array(f["PartType1/Velocities"])
e_kin_particles = 0.5*v**2
e_kin_total = np.sum(e_kin_particles)
e_kin_array = np.append(e_kin_array,e_kin_total)
e_total_array = np.append(e_total_array,e_kin_total+e_pot_total)
print snap_numbers.shape
print e_kin_array.shape
plt.plot(e_kin_array,label = "Kinetic Energy")
plt.plot(e_pot_array,label = "Potential Energy")
plt.plot(e_total_array,label = "Total Energy")
plt.legend(loc = "upper right")
plot_filename = "%s/energy_plot.png" %plot_dir
plt.savefig(plot_filename,format = "png")
plt.close()
import h5py as h5
import numpy as np
import matplotlib.pyplot as plt
n_snaps = 102
directory = '../'
snap_name_base = "Isothermal_"
plot_directory = "random_orbits_3d"
plot_filename_base = "snap"
for i in range(n_snaps):
snap_number = "%03d" %i
filename = directory + snap_name_base + snap_number + ".hdf5"
f = h5.File(filename)
#find the box size
header = f["Header"]
box_size = header.attrs["BoxSize"][0]
pos = np.array(f["PartType1/Coordinates"])
x = pos[:,0]
y = pos[:,1]
z = pos[:,2]
#x_projection
plt.plot(y,z,'bo')
plt.xlabel("y")
plt.ylabel("z")
plt.xlim((0.0,box_size))
plt.ylim((0.0,box_size))
plot_dir = "%s/projection_x" %plot_directory
plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
plt.savefig(plot_filename,format = "png")
plt.close()
#y_projection
plt.plot(x,z,'bo')
plt.xlabel("x")
plt.ylabel("z")
plt.xlim((0.0,box_size))
plt.ylim((0.0,box_size))
plot_dir = "%s/projection_y" %plot_directory
plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
plt.savefig(plot_filename,format = "png")
plt.close()
#x_projection
plt.plot(x,y,'bo')
plt.xlabel("x")
plt.ylabel("y")
plt.xlim((0.0,box_size))
plt.ylim((0.0,box_size))
plot_dir = "%s/projection_z" %plot_directory
plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
plt.savefig(plot_filename,format = "png")
plt.close()
......@@ -4,7 +4,7 @@
if [ ! -e Isothermal.hdf5 ]
then
echo "Generating initial conditions for the isothermal potential box example..."
python makeIC.py 1000 1
python makeIC.py 10000 1 1
fi
../../swift -g -t 2 isothermal.yml
import numpy as np
import matplotlib.pyplot as plt
import h5py as h5
n_snaps = 102
directory = '../'
snap_name_base = "Isothermal_"
plot_dir = "random_orbits_1d"
plot_filename_base = "snap"
v_rot = 200.
r_0 = 100.
e_kin_array = []
e_pot_array = []
e_total_array = []
snap_numbers = np.linspace(0,n_snaps)
for i in range(n_snaps):
snap_number = "%03d" %i
filename = directory + snap_name_base + snap_number + ".hdf5"
f = h5.File(filename)
header = f["Header"]
n_parts = header.attrs["NumPart_ThisFile"][1]
pos = np.array(f["PartType1/Coordinates"])
x = pos[:,0]
y = pos[:,1]
z = pos[:,2]
r = np.sqrt(x**2 + y**2 + z**2)
e_pot_particles = v_rot**2 * np.log(r/r_0)
e_pot_total = np.sum(e_pot_particles)
e_pot_array = np.append(e_pot_array,e_pot_total)
v = np.array(f["PartType1/Velocities"])
e_kin_particles = 0.5*v**2
e_kin_total = np.sum(e_kin_particles)
e_kin_array = np.append(e_kin_array,e_kin_total)
e_total_array = np.append(e_total_array,e_kin_total+e_pot_total)
print snap_numbers.shape
print e_kin_array.shape
plt.plot(e_kin_array,label = "Kinetic Energy")
plt.plot(e_pot_array,label = "Potential Energy")
plt.plot(e_total_array,label = "Total Energy")
plt.legend(loc = "upper right")
plot_filename = "%s/energy_plot.png" %plot_dir
plt.savefig(plot_filename,format = "png")
plt.close()
import h5py as h5
import numpy as np
import matplotlib.pyplot as plt
n_snaps = 100
directory = '../'
snap_name_base = "uniformBox_"
plot_directory = "./"
plot_filename_base = "snap"
for i in range(n_snaps):
snap_number = "%03d" %i
filename = directory + snap_name_base + snap_number + ".hdf5"
f = h5.File(filename)
#find the box size
header = f["Header"]
box_size = header.attrs["BoxSize"][0]
pos = np.array(f["PartType0/Coordinates"])
x = pos[:,0]
y = pos[:,1]
z = pos[:,2]
#x_projection
plt.plot(y,z,'bo')
plt.xlabel("y")
plt.ylabel("z")
plt.xlim((0.0,box_size))
plt.ylim((0.0,box_size))
plot_dir = "%s/projection_x" %plot_directory
plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
plt.savefig(plot_filename,format = "png")
plt.close()
#y_projection
plt.plot(x,z,'bo')
plt.xlabel("x")
plt.ylabel("z")
plt.xlim((0.0,box_size))
plt.ylim((0.0,box_size))
plot_dir = "%s/projection_y" %plot_directory
plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
plt.savefig(plot_filename,format = "png")
plt.close()
#x_projection
plt.plot(x,y,'bo')
plt.xlabel("x")
plt.ylabel("y")
plt.xlim((0.0,box_size))
plt.ylim((0.0,box_size))
plot_dir = "%s/projection_z" %plot_directory
plot_filename = "%s/%s_%03d.png" %(plot_dir,plot_filename_base,i)
plt.savefig(plot_filename,format = "png")
plt.close()
......@@ -7,4 +7,4 @@ then
python makeIC.py 100
fi
../swift -s -t 16 uniformBox.yml
../swift -s -C -t 16 uniformBox.yml
......@@ -58,6 +58,7 @@ void print_help_message() {
printf("Valid options are:\n");
printf(" %2s %8s %s\n", "-a", "", "Pin runners using processor affinity");
printf(" %2s %8s %s\n", "-c", "", "Run with cosmological time integration");
printf(" %2s %8s %s\n", "-C", "", "Run with cooling");
printf(
" %2s %8s %s\n", "-d", "",
"Dry run. Read the parameter file, allocate memory but does not read ");
......@@ -154,7 +155,7 @@ int main(int argc, char *argv[]) {
/* Parse the parameters */
int c;
while ((c = getopt(argc, argv, "acdef:gGhn:st:v:y:")) != -1) switch (c) {
while ((c = getopt(argc, argv, "acCdef:gGhn:st:v:y:")) != -1) switch (c) {
case 'a':
with_aff = 1;
break;
......
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