Commit 32298b94 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'master' into autotools-update

parents 27174568 907b0ce7
...@@ -22,8 +22,6 @@ doc/Doxyfile ...@@ -22,8 +22,6 @@ doc/Doxyfile
examples/swift examples/swift
examples/swift_mpi examples/swift_mpi
examples/swift_fixdt
examples/swift_fixdt_mpi
examples/*.xmf examples/*.xmf
examples/used_parameters.yml examples/used_parameters.yml
examples/energy.txt examples/energy.txt
...@@ -80,9 +78,12 @@ tests/testRiemannHLLC ...@@ -80,9 +78,12 @@ tests/testRiemannHLLC
tests/testMatrixInversion tests/testMatrixInversion
theory/latex/swift.pdf theory/latex/swift.pdf
theory/kernel/kernels.pdf theory/SPH/Kernels/kernels.pdf
theory/kernel/kernel_derivatives.pdf theory/SPH/Kernels/kernel_derivatives.pdf
theory/kernel/kernel_definitions.pdf theory/SPH/Kernels/kernel_definitions.pdf
theory/SPH/Flavours/sph_flavours.pdf
theory/SPH/EoS/eos.pdf
theory/SPH/*.pdf
theory/paper_pasc/pasc_paper.pdf theory/paper_pasc/pasc_paper.pdf
m4/libtool.m4 m4/libtool.m4
......
...@@ -83,39 +83,65 @@ SWIFT depends on a number of third party libraries that should be available ...@@ -83,39 +83,65 @@ SWIFT depends on a number of third party libraries that should be available
before you can build it. before you can build it.
HDF5: a HDF5 library (v. 1.8.x or higher) is required to read and write - HDF5: a HDF5 library (v. 1.8.x or higher) is required to read and
particle data. One of the commands "h5cc" or "h5pcc" should be write particle data. One of the commands "h5cc" or "h5pcc"
available. If "h5pcc" is located them a parallel HDF5 built for the version should be available. If "h5pcc" is located them a parallel
of MPI located should be provided. If the command is not available then it HDF5 built for the version of MPI located should be
can be located using the "--with-hfd5" configure option. The value should provided. If the command is not available then it can be
be the full path to the "h5cc" or "h5pcc" commands. located using the "--with-hfd5" configure option. The value
should be the full path to the "h5cc" or "h5pcc" commands.
MPI: an optional MPI library that fully supports MPI_THREAD_MULTIPLE. - MPI: to run on more than one node an MPI library that fully
Before running configure the "mpirun" command should be available in the supports MPI_THREAD_MULTIPLE. Before running configure the
shell. If your command isn't called "mpirun" then define the "MPIRUN" "mpirun" command should be available in the shell. If your
environment variable, either in the shell or when running configure. command isn't called "mpirun" then define the "MPIRUN"
environment variable, either in the shell or when running
configure.
The MPI compiler can be controlled using the MPICC variable, much like The MPI compiler can be controlled using the MPICC variable,
the CC one. Use this when your MPI compiler has a none-standard name. much like the CC one. Use this when your MPI compiler has a
none-standard name.
METIS: a build of the METIS library can be optionally used to optimize the - libtool: The build system relies on libtool.
load between MPI nodes (requires an MPI library). This should be found in
the standard installation directories, or pointed at using the
"--with-metis" configuration option. In this case the top-level
installation directory of the METIS build should be given. Note to use
METIS you should at least supply "--with-metis".
libNUMA: a build of the NUMA library can be used to pin the threads to Optional Dependencies
the physical core of the machine SWIFT is running on. This is not always =====================
necessary as the OS scheduler may do a good job at distributing the threads
among the different cores on each computing node.
DOXYGEN: the doxygen library is required to create the SWIFT API - METIS: a build of the METIS library can be optionally used to
documentation. optimize the load between MPI nodes (requires an MPI
library). This should be found in the standard installation
directories, or pointed at using the "--with-metis"
configuration option. In this case the top-level
installation directory of the METIS build should be
given. Note to use METIS you should at least supply
"--with-metis".
- libNUMA: a build of the NUMA library can be used to pin the threads
to the physical core of the machine SWIFT is running
on. This is not always necessary as the OS scheduler may
do a good job at distributing the threads among the
different cores on each computing node.
- TCMalloc: a build of the TCMalloc library (part of gperftools) can
be used to obtain faster allocations than the standard C
malloc function part of glibc. The option "-with-tcmalloc"
should be passed to the configuration script to use it.
- gperftools: a build of gperftools can be used to obtain good
profiling of the code. The option "-with-profiler"
needs to be passed to the configuration script to use
it.
- DOXYGEN: the doxygen library is required to create the SWIFT API
documentation.
......
...@@ -13,8 +13,6 @@ See INSTALL.swift for install instructions. ...@@ -13,8 +13,6 @@ See INSTALL.swift for install instructions.
Usage: swift [OPTION]... PARAMFILE Usage: swift [OPTION]... PARAMFILE
swift_mpi [OPTION]... PARAMFILE swift_mpi [OPTION]... PARAMFILE
swift_fixdt [OPTION]... PARAMFILE
swift_fixdt_mpi [OPTION]... PARAMFILE
Valid options are: Valid options are:
-a Pin runners using processor affinity -a Pin runners using processor affinity
......
...@@ -170,6 +170,18 @@ if test "x$enable_debug" = "xyes"; then ...@@ -170,6 +170,18 @@ if test "x$enable_debug" = "xyes"; then
fi fi
fi fi
# Check if task debugging is on.
AC_ARG_ENABLE([task-debugging],
[AS_HELP_STRING([--enable-task-debugging],
[Store task timing information and generate task dump files @<:@yes/no@:>@]
)],
[enable_task_debugging="$enableval"],
[enable_task_debugging="no"]
)
if test "$enable_task_debugging" = "yes"; then
AC_DEFINE([SWIFT_DEBUG_TASKS],1,[Enable task debugging])
fi
# Define HAVE_POSIX_MEMALIGN if it works. # Define HAVE_POSIX_MEMALIGN if it works.
AX_FUNC_POSIX_MEMALIGN AX_FUNC_POSIX_MEMALIGN
...@@ -355,7 +367,7 @@ fi ...@@ -355,7 +367,7 @@ fi
AC_SUBST([TCMALLOC_LIBS]) AC_SUBST([TCMALLOC_LIBS])
AM_CONDITIONAL([HAVETCMALLOC],[test -n "$TCMALLOC_LIBS"]) AM_CONDITIONAL([HAVETCMALLOC],[test -n "$TCMALLOC_LIBS"])
# Check for -lprofiler usually part of the gpreftools along with tcmalloc. # Check for -lprofiler usually part of the gperftools along with tcmalloc.
have_profiler="no" have_profiler="no"
AC_ARG_WITH([profiler], AC_ARG_WITH([profiler],
[AS_HELP_STRING([--with-profiler], [AS_HELP_STRING([--with-profiler],
...@@ -595,6 +607,7 @@ AC_MSG_RESULT([ ...@@ -595,6 +607,7 @@ AC_MSG_RESULT([
CPU profiler : $have_profiler CPU profiler : $have_profiler
Hydro scheme : $with_hydro Hydro scheme : $with_hydro
Dimensionality : $with_dimension Dimensionality : $with_dimension
Task debugging : $enable_task_debugging
]) ])
# Generate output. # Generate output.
......
# Define the system of units to use internally. # Define the system of units to use internally.
InternalUnitSystem: InternalUnitSystem:
UnitMass_in_cgs: 2.0e33 # Solar masses UnitMass_in_cgs: 2.0e33 # Solar masses
UnitLength_in_cgs: 3.01e21 # Kilparsecs UnitLength_in_cgs: 3.0857e21 # Kiloparsecs
UnitVelocity_in_cgs: 1.0e5 # Time unit is cooling time UnitVelocity_in_cgs: 1.0e5 # Time unit is cooling time
UnitCurrent_in_cgs: 1 # Amperes UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin UnitTemp_in_cgs: 1 # Kelvin
...@@ -9,9 +9,9 @@ InternalUnitSystem: ...@@ -9,9 +9,9 @@ InternalUnitSystem:
# Parameters governing the time integration # Parameters governing the time integration
TimeIntegration: TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units). time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 1.0 # The end time of the simulation (in internal units). time_end: 4. # The end time of the simulation (in internal units).
dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units). dt_min: 1e-4 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units). dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots # Parameters governing the snapshots
Snapshots: Snapshots:
...@@ -27,7 +27,6 @@ Statistics: ...@@ -27,7 +27,6 @@ Statistics:
SPH: SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel). resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours. delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
max_smoothing_length: 0.1 # Maximal smoothing length allowed (in internal units).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration. CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
# Parameters related to the initial conditions # Parameters related to the initial conditions
...@@ -36,9 +35,8 @@ InitialConditions: ...@@ -36,9 +35,8 @@ InitialConditions:
# Dimensionless pre-factor for the time-step condition # Dimensionless pre-factor for the time-step condition
LambdaCooling: LambdaCooling:
lambda: 0.0 # Cooling rate (in cgs units) lambda_cgs: 1.0e-22 # Cooling rate (in cgs units)
minimum_temperature: 1.0e4 # Minimal temperature (Kelvin) minimum_temperature: 1.0e4 # Minimal temperature (Kelvin)
mean_molecular_weight: 0.59 # Mean molecular weight mean_molecular_weight: 0.59 # Mean molecular weight
hydrogen_mass_abundance: 0.75 # Hydrogen mass abundance (dimensionless) hydrogen_mass_abundance: 0.75 # Hydrogen mass abundance (dimensionless)
cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
...@@ -13,7 +13,7 @@ m_p = 1.67e-24 #proton mass ...@@ -13,7 +13,7 @@ m_p = 1.67e-24 #proton mass
#initial conditions set in makeIC.py #initial conditions set in makeIC.py
rho = 3.2e3 rho = 3.2e3
P = 4.5e6 P = 4.5e6
n_H_cgs = 0.0001 #n_H_cgs = 0.0001
gamma = 5./3. gamma = 5./3.
T_init = 1.0e5 T_init = 1.0e5
...@@ -24,7 +24,7 @@ unit_mass = units.attrs["Unit mass in cgs (U_M)"] ...@@ -24,7 +24,7 @@ unit_mass = units.attrs["Unit mass in cgs (U_M)"]
unit_length = units.attrs["Unit length in cgs (U_L)"] unit_length = units.attrs["Unit length in cgs (U_L)"]
unit_time = units.attrs["Unit time in cgs (U_t)"] unit_time = units.attrs["Unit time in cgs (U_t)"]
parameters = f["Parameters"] parameters = f["Parameters"]
cooling_lambda = float(parameters.attrs["LambdaCooling:lambda"]) cooling_lambda = float(parameters.attrs["LambdaCooling:lambda_cgs"])
min_T = float(parameters.attrs["LambdaCooling:minimum_temperature"]) min_T = float(parameters.attrs["LambdaCooling:minimum_temperature"])
mu = float(parameters.attrs["LambdaCooling:mean_molecular_weight"]) mu = float(parameters.attrs["LambdaCooling:mean_molecular_weight"])
X_H = float(parameters.attrs["LambdaCooling:hydrogen_mass_abundance"]) X_H = float(parameters.attrs["LambdaCooling:hydrogen_mass_abundance"])
...@@ -32,50 +32,65 @@ X_H = float(parameters.attrs["LambdaCooling:hydrogen_mass_abundance"]) ...@@ -32,50 +32,65 @@ X_H = float(parameters.attrs["LambdaCooling:hydrogen_mass_abundance"])
#get number of particles #get number of particles
header = f["Header"] header = f["Header"]
n_particles = header.attrs["NumPart_ThisFile"][0] n_particles = header.attrs["NumPart_ThisFile"][0]
#read energy and time arrays #read energy and time arrays
array = np.genfromtxt(stats_filename,skip_header = 1) array = np.genfromtxt(stats_filename,skip_header = 1)
time = array[:,0] time = array[:,0]
total_energy = array[:,2] kin_plus_therm = array[:,2]
radiated = array[:,6]
total_mass = array[:,1] total_mass = array[:,1]
#ignore first row where there are just zeros
time = time[1:] time = time[1:]
total_energy = total_energy[1:] kin_plus_therm = kin_plus_therm[1:]
radiated = radiated[1:]
total_mass = total_mass[1:] total_mass = total_mass[1:]
total_energy = kin_plus_therm + radiated
initial_energy = total_energy[0]
#conversions to cgs #conversions to cgs
rho_cgs = rho * unit_mass / (unit_length)**3 rho_cgs = rho * unit_mass / (unit_length)**3
time_cgs = time * unit_time time_cgs = time * unit_time
u_init_cgs = total_energy[0]/(total_mass[0]) * unit_length**2 / (unit_time)**2 initial_energy_cgs = initial_energy/total_mass[0] * unit_length**2 / (unit_time)**2
n_H_cgs = X_H * rho_cgs / m_p
#find the energy floor #find the energy floor
print min_T
u_floor_cgs = k_b * min_T / (mu * m_p * (gamma - 1.)) u_floor_cgs = k_b * min_T / (mu * m_p * (gamma - 1.))
#find analytic solution #find analytic solution
analytic_time = np.linspace(time_cgs[0],time_cgs[-1],1000) analytic_time_cgs = np.linspace(0,time_cgs[-1],1000)
print time_cgs[1]
print analytic_time[1]
du_dt_cgs = -cooling_lambda * n_H_cgs**2 / rho_cgs du_dt_cgs = -cooling_lambda * n_H_cgs**2 / rho_cgs
u_analytic = du_dt_cgs*(analytic_time - analytic_time[0]) + u_init_cgs u_analytic_cgs = du_dt_cgs*analytic_time_cgs + initial_energy_cgs
cooling_time = u_init_cgs/(-du_dt_cgs) cooling_time_cgs = initial_energy_cgs/(-du_dt_cgs)
#rescale energy to initial energy
total_energy /= total_energy[0] for i in range(u_analytic_cgs.size):
u_analytic /= u_init_cgs if u_analytic_cgs[i]<u_floor_cgs:
u_floor_cgs /= u_init_cgs u_analytic_cgs[i] = u_floor_cgs
# plot_title = r"$\Lambda \, = \, %1.1g \mathrm{erg}\mathrm{cm^3}\mathrm{s^{-1}} \, \, T_{init} = %1.1g\mathrm{K} \, \, T_{floor} = %1.1g\mathrm{K} \, \, n_H = %1.1g\mathrm{cm^{-3}}$" %(cooling_lambda,T_init,T_floor,n_H)
# plot_filename = "energy_plot_creasey_no_cooling_T_init_1p0e5_n_H_0p1.png" #rescale analytic solution
#analytic_solution = np.zeros(n_snaps-1) u_analytic = u_analytic_cgs/initial_energy_cgs
for i in range(u_analytic.size):
if u_analytic[i]<u_floor_cgs: #put time in units of cooling_time
u_analytic[i] = u_floor_cgs time=time_cgs/cooling_time_cgs
plt.plot(time_cgs,total_energy,'k',label = "Numerical solution") analytic_time = analytic_time_cgs/cooling_time_cgs
plt.plot(analytic_time,u_analytic,'--r',lw = 2.0,label = "Analytic Solution")
plt.plot((cooling_time,cooling_time),(0,1),'b',label = "Cooling time") #rescale (numerical) energy by initial energy
plt.plot((time_cgs[0],time_cgs[0]),(0,1),'m',label = "First output") radiated /= initial_energy
plt.title(r"$n_H = %1.1e \, \mathrm{cm}^{-3}$" %n_H_cgs) kin_plus_therm /= initial_energy
plt.xlabel("Time (seconds)") total_energy = kin_plus_therm + radiated
plt.ylabel("Energy/Initial energy") plt.plot(time,kin_plus_therm,'kd',label = "Kinetic + thermal energy")
plt.plot(time,radiated,'bo',label = "Radiated energy")
plt.plot(time,total_energy,'g',label = "Total energy")
plt.plot(analytic_time,u_analytic,'r',lw = 2.0,label = "Analytic Solution")
#plt.plot(analytic_time,1-u_analytic,'k',lw = 2.0)
#plt.plot((cooling_time,cooling_time),(0,1),'b',label = "Cooling time")
#plt.plot((time[1]-time_cgs[0],time_cgs[1]-time_cgs[0]),(0,1),'m',label = "First output")
#plt.title(r"$n_H = %1.1e \, \mathrm{cm}^{-3}$" %n_H_cgs)
plt.xlabel("Time / cooling time")
plt.ylabel("Energy / Initial energy")
#plt.ylim(0,1.1)
plt.ylim(0.999,1.001) plt.ylim(0.999,1.001)
plt.xlim(0,min(10.0*cooling_time,time_cgs[-1])) #plt.xlim(0,min(10,time[-1]))
plt.legend(loc = "upper right") plt.legend(loc = "upper right")
if (int(sys.argv[1])==0): if (int(sys.argv[1])==0):
plt.show() plt.show()
......
...@@ -27,10 +27,10 @@ from numpy import * ...@@ -27,10 +27,10 @@ from numpy import *
# Parameters # Parameters
periodic= 1 # 1 For periodic box periodic= 1 # 1 For periodic box
boxSize = 1 #1 kiloparsec boxSize = 1 # 1 kiloparsec
L = int(sys.argv[1]) # Number of particles along one axis L = int(sys.argv[1]) # Number of particles along one axis
rho = 3.2e3 # Density in code units (0.01 hydrogen atoms per cm^3) rho = 3.2e3 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
P = 4.5e6 # Pressure in code units (at 10^5K) P = 4.5e6 # Pressure in code units (at 10^5K)
gamma = 5./3. # Gas adiabatic index gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "coolingBox.hdf5" fileName = "coolingBox.hdf5"
...@@ -63,9 +63,9 @@ grp.attrs["PeriodicBoundariesOn"] = periodic ...@@ -63,9 +63,9 @@ grp.attrs["PeriodicBoundariesOn"] = periodic
#Units #Units
grp = file.create_group("/Units") grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = 3.08e21 grp.attrs["Unit length in cgs (U_L)"] = 3.0857e21
grp.attrs["Unit mass in cgs (U_M)"] = 2.0e33 grp.attrs["Unit mass in cgs (U_M)"] = 2.0e33
grp.attrs["Unit time in cgs (U_t)"] = 3.08e16 grp.attrs["Unit time in cgs (U_t)"] = 3.0857e16
grp.attrs["Unit current in cgs (U_I)"] = 1. grp.attrs["Unit current in cgs (U_I)"] = 1.
grp.attrs["Unit temperature in cgs (U_T)"] = 1. grp.attrs["Unit temperature in cgs (U_T)"] = 1.
......
...@@ -5,6 +5,10 @@ echo "Generating initial conditions for the cooling box example..." ...@@ -5,6 +5,10 @@ echo "Generating initial conditions for the cooling box example..."
python makeIC.py 10 python makeIC.py 10
../swift -s -t 1 coolingBox.yml -C 2>&1 | tee output.log ../swift -s -C -t 16 coolingBox.yml
#-C 2>&1 | tee output.log
python energy_plot.py 0 python energy_plot.py 0
#python test_energy_conservation.py 0
import numpy as np
import matplotlib.pyplot as plt
import h5py as h5
import sys
stats_filename = "./energy.txt"
snap_filename = "coolingBox_000.hdf5"
#plot_dir = "./"
n_snaps = 41
time_end = 4.0
dt_snap = 0.1
#some constants in cgs units
k_b = 1.38E-16 #boltzmann
m_p = 1.67e-24 #proton mass
#initial conditions set in makeIC.py
rho = 4.8e3
P = 4.5e6
#n_H_cgs = 0.0001
gamma = 5./3.
T_init = 1.0e5
#find the sound speed
#Read the units parameters from the snapshot
f = h5.File(snap_filename,'r')
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)"]
parameters = f["Parameters"]
cooling_lambda = float(parameters.attrs["LambdaCooling:lambda_cgs"])
min_T = float(parameters.attrs["LambdaCooling:minimum_temperature"])
mu = float(parameters.attrs["LambdaCooling:mean_molecular_weight"])
X_H = float(parameters.attrs["LambdaCooling:hydrogen_mass_abundance"])
#get number of particles
header = f["Header"]
n_particles = header.attrs["NumPart_ThisFile"][0]
#read energy and time arrays
array = np.genfromtxt(stats_filename,skip_header = 1)
time = array[:,0]
total_energy = array[:,2]
total_mass = array[:,1]
time = time[1:]
total_energy = total_energy[1:]
total_mass = total_mass[1:]
#conversions to cgs
rho_cgs = rho * unit_mass / (unit_length)**3
time_cgs = time * unit_time
u_init_cgs = total_energy[0]/(total_mass[0]) * unit_length**2 / (unit_time)**2
n_H_cgs = X_H * rho_cgs / m_p
#find the sound speed in cgs
c_s = np.sqrt((gamma - 1.)*k_b*T_init/(mu*m_p))
#assume box size is unit length
sound_crossing_time = unit_length/c_s
print "Sound speed = %g cm/s" %c_s
print "Sound crossing time = %g s" %sound_crossing_time
#find the energy floor
u_floor_cgs = k_b * min_T / (mu * m_p * (gamma - 1.))
#find analytic solution
analytic_time_cgs = np.linspace(time_cgs[0],time_cgs[-1],1000)
du_dt_cgs = -cooling_lambda * n_H_cgs**2 / rho_cgs
u_analytic = du_dt_cgs*(analytic_time_cgs - analytic_time_cgs[0]) + u_init_cgs
cooling_time = u_init_cgs/(-du_dt_cgs)
#put time in units of sound crossing time
time=time_cgs/sound_crossing_time
analytic_time = analytic_time_cgs/sound_crossing_time
#rescale energy to initial energy
total_energy /= total_energy[0]
u_analytic /= u_init_cgs
u_floor_cgs /= u_init_cgs
# plot_title = r"$\Lambda \, = \, %1.1g \mathrm{erg}\mathrm{cm^3}\mathrm{s^{-1}} \, \, T_{init} = %1.1g\mathrm{K} \, \, T_{floor} = %1.1g\mathrm{K} \, \, n_H = %1.1g\mathrm{cm^{-3}}$" %(cooling_lambda,T_init,T_floor,n_H)
# plot_filename = "energy_plot_creasey_no_cooling_T_init_1p0e5_n_H_0p1.png"
#analytic_solution = np.zeros(n_snaps-1)
for i in range(u_analytic.size):
if u_analytic[i]<u_floor_cgs:
u_analytic[i] = u_floor_cgs
plt.plot(time-time[0],total_energy,'k',label = "Numerical solution from energy.txt")
plt.plot(analytic_time-analytic_time[0],u_analytic,'r',lw = 2.0,label = "Analytic Solution")
#now get energies from the snapshots
snapshot_time = np.linspace(0,time_end,num = n_snaps)
snapshot_time = snapshot_time[1:]
snapshot_time_cgs = snapshot_time * unit_time
snapshot_time = snapshot_time_cgs/ sound_crossing_time
snapshot_time -= snapshot_time[0]
snapshot_energy = np.zeros(n_snaps)
for i in range(0,n_snaps):
snap_filename = "coolingBox_%03d.hdf5" %i
f = h5.File(snap_filename,'r')
snapshot_internal_energy_array = np.array(f["PartType0/InternalEnergy"])
total_internal_energy = np.sum(snapshot_internal_energy_array)
velocity_array = np.array(f["PartType0/Velocities"])
total_kinetic_energy = 0.5*np.sum(velocity_array**2)
snapshot_energy[i] = total_internal_energy + total_kinetic_energy
snapshot_energy/=snapshot_energy[0]
snapshot_energy = snapshot_energy[1:]
plt.plot(snapshot_time,snapshot_energy,'bd',label = "Numerical solution from snapshots")
#plt.title(r"$n_H = %1.1e \, \mathrm{cm}^{-3}$" %n_H_cgs)
plt.xlabel("Time (sound crossing time)")
plt.ylabel("Energy/Initial energy")
plt.ylim(0.99,1.01)
#plt.xlim(0,min(10,time[-1]))
plt.legend(loc = "upper right")
if (int(sys.argv[1])==0):
plt.show()
else:
plt.savefig(full_plot_filename,format = "png")
plt.close()
To make the initial conditions we distribute gas particles randomly in
a cube with a side length twice that of the virial radius. The density
profile of the gas is proportional to r^(-2) where r is the distance
from the centre of the cube.
The parameter v_rot (in makeIC.py and cooling.yml) sets the circular
velocity of the halo, and by extension, the viral radius, viral mass,
and the internal energy of the gas such that hydrostatic equilibrium
is achieved.
While the system is initially in hydrostatic equilibrium, the cooling
of the gas means that the halo will collapse.
To run this example, make such that the code is compiled with either
the isothermal potential or softened isothermal potential, and
'const_lambda' cooling, set in src/const.h. In the latter case, a
(small) value of epsilon needs to be set in cooling.yml. 0.1 kpc
should work well.
The plotting scripts produce a plot of the density, internal energy
and radial velocity profile for each
snapshot. test_energy_conservation.py shows the evolution of energy
with time. These can be used to check if the example has run properly.
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.9885e39 # 10^6 solar masses
UnitLength_in_cgs: 3.0856776e21 # Kiloparsecs
UnitVelocity_in_cgs: 1e5 # Kilometres per second
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: 10.0 # The end time of the simulation (in internal units).
dt_min: 1e-5 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
# Parameters governing the snapshots