...
 
Commits (1)
This example generates a set of particles in an isothermal potential
and follows their orbits. IDL scripts verify the conservation of
energy and angular momentum.
import matplotlib
matplotlib.use("Agg")
from pylab import *
import h5py
# Plot parameters
params = {'axes.labelsize': 10,
'axes.titlesize': 10,
'font.size': 12,
'legend.fontsize': 12,
'xtick.labelsize': 10,
'ytick.labelsize': 10,
'text.usetex': True,
'figure.figsize' : (3.15,3.15),
'figure.subplot.left' : 0.145,
'figure.subplot.right' : 0.99,
'figure.subplot.bottom' : 0.11,
'figure.subplot.top' : 0.99,
'figure.subplot.wspace' : 0.15,
'figure.subplot.hspace' : 0.12,
'lines.markersize' : 6,
'lines.linewidth' : 3.,
'text.latex.unicode': True
}
rcParams.update(params)
rc('font',**{'family':'sans-serif','sans-serif':['Times']})
import numpy as np
import h5py as h5
import sys
# File containing the total energy
stats_filename = "./energy.txt"
# First snapshot
snap_filename = "Isothermal_0000.hdf5"
f = h5.File(snap_filename,'r')
# Read the units parameters from the snapshot
units = f["InternalCodeUnits"]
unit_mass = units.attrs["Unit mass in cgs (U_M)"]
unit_length = units.attrs["Unit length in cgs (U_L)"]
unit_time = units.attrs["Unit time in cgs (U_t)"]
# Read the header
header = f["Header"]
box_size = float(header.attrs["BoxSize"][0])
# Read the properties of the potential
parameters = f["Parameters"]
R200 = 100
Vrot = float(parameters.attrs["IsothermalPotential:vrot"])
centre = [box_size/2, box_size/2, box_size/2]
f.close()
# Read the statistics summary
file_energy = np.loadtxt("energy.txt")
time_stats = file_energy[:,0]
E_kin_stats = file_energy[:,3]
E_pot_stats = file_energy[:,5]
E_tot_stats = E_kin_stats + E_pot_stats
# Read the snapshots
time_snap = np.zeros(402)
E_kin_snap = np.zeros(402)
E_pot_snap = np.zeros(402)
E_tot_snap = np.zeros(402)
Lz_snap = np.zeros(402)
# Read all the particles from the snapshots
for i in range(402):
snap_filename = "Isothermal_%0.4d.hdf5"%i
f = h5.File(snap_filename,'r')
pos_x = f["PartType3/Coordinates"][:,0]
pos_y = f["PartType3/Coordinates"][:,1]
pos_z = f["PartType3/Coordinates"][:,2]
vel_x = f["PartType3/Velocities"][:,0]
vel_y = f["PartType3/Velocities"][:,1]
vel_z = f["PartType3/Velocities"][:,2]
mass = f["/PartType3/Masses"][:]
r = np.sqrt((pos_x[:] - centre[0])**2 + (pos_y[:] - centre[1])**2 + (pos_z[:] - centre[2])**2)
Lz = (pos_x[:] - centre[0]) * vel_y[:] - (pos_y[:] - centre[1]) * vel_x[:]
time_snap[i] = f["Header"].attrs["Time"]
E_kin_snap[i] = np.sum(0.5 * mass * (vel_x[:]**2 + vel_y[:]**2 + vel_z[:]**2))
E_pot_snap[i] = np.sum(mass * Vrot**2 * log(r))
E_tot_snap[i] = E_kin_snap[i] + E_pot_snap[i]
Lz_snap[i] = np.sum(Lz)
# Plot energy evolution
figure()
plot(time_stats, E_kin_stats, "r-", lw=0.5, label="Kinetic energy")
plot(time_stats, E_pot_stats, "g-", lw=0.5, label="Potential energy")
plot(time_stats, E_tot_stats, "k-", lw=0.5, label="Total energy")
plot(time_snap[::10], E_kin_snap[::10], "rD", lw=0.5, ms=2)
plot(time_snap[::10], E_pot_snap[::10], "gD", lw=0.5, ms=2)
plot(time_snap[::10], E_tot_snap[::10], "kD", lw=0.5, ms=2)
legend(loc="center right", fontsize=8, frameon=False, handlelength=3, ncol=1)
xlabel("${\\rm{Time}}$", labelpad=0)
ylabel("${\\rm{Energy}}$",labelpad=0)
xlim(0, 8)
savefig("energy.png", dpi=200)
# Plot angular momentum evolution
figure()
plot(time_snap, Lz_snap, "k-", lw=0.5, ms=2)
xlabel("${\\rm{Time}}$", labelpad=0)
ylabel("${\\rm{Angular~momentum}}$",labelpad=0)
xlim(0, 8)
savefig("angular_momentum.png", dpi=200)
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e33 # M_sun
UnitLength_in_cgs: 3.08567758e21 # kpc
UnitVelocity_in_cgs: 1e5 # km/s
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# 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).
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).
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-3 # Time between statistics output
# Parameters governing the snapshots
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)
# Parameters related to the initial conditions
InitialConditions:
file_name: Isothermal.hdf5 # The file to read
periodic: 1
shift: [200.,200.,200.] # Shift all particles to be in the potential
# External potential parameters
IsothermalPotential:
useabspos: 0 # Whether to use absolute position (1) or relative potential to centre of box (0)
position: [0.,0.,0.]
vrot: 200. # rotation speed of isothermal potential in internal units
timestep_mult: 0.01 # controls time step
epsilon: 0. # No softening at the centre of the halo
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##############################################################################
import h5py
import sys
import numpy
import math
# 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]
# all particles move in the xy plane, and start at y=0
# physical constants in cgs
NEWTON_GRAVITY_CGS = 6.67408e-8
SOLAR_MASS_IN_CGS = 1.98848e33
PARSEC_IN_CGS = 3.08567758e18
YEAR_IN_CGS = 3.15569252e7
# choice of units
const_unit_length_in_cgs = (1000*PARSEC_IN_CGS)
const_unit_mass_in_cgs = (SOLAR_MASS_IN_CGS)
const_unit_velocity_in_cgs = (1e5)
print("UnitMass_in_cgs: ", const_unit_mass_in_cgs)
print("UnitLength_in_cgs: ", const_unit_length_in_cgs)
print("UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs)
# rotation speed of isothermal potential [km/s]
vrot_kms = 200.
# derived units
const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
const_G = ((NEWTON_GRAVITY_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
print('G=', const_G)
vrot = vrot_kms * 1e5 / const_unit_velocity_in_cgs
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 400. # [kpc]
Radius = 100. # maximum radius of particles [kpc]
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[
L = N**(1./3.)
fileName = "Isothermal.hdf5"
#---------------------------------------------------
numPart = N
mass = 1
#--------------------------------------------------
#File
file = h5py.File(fileName, 'w')
#Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs
grp.attrs["Unit time in cgs (U_t)"] = const_unit_length_in_cgs / const_unit_velocity_in_cgs
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"] = [0, 0, 0, numPart, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [0, 0, 0, numPart, 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
# set seed for random number
numpy.random.seed(1234)
#Particle group
grp1 = file.create_group("/PartType3")
#generate particle positions
radius = Radius * (numpy.random.rand(N))**(1./3.)
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
#
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]
ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
v = numpy.zeros(1)
m = numpy.full((numPart, ), mass, dtype='f')
ds = grp1.create_dataset('Masses', (numPart,), 'f')
ds[()] = m
m = numpy.zeros(1)
ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False)
ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
ds[()] = ids
ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
ds[()] = r
file.close()
#!/bin/bash
# Generate the initial conditions if they are not present.
if [ ! -e Isothermal.hdf5 ]
then
echo "Generating initial conditions for the isothermal potential box example..."
python makeIC.py 1000 0
fi
rm -rf Isothermal_*.hdf5
../../swift --sinks --external-gravity --threads=1 isothermal.yml 2>&1 | tee output.log
python energy_plot.py