Commit 3278ea21 authored by James Willis's avatar James Willis
Browse files

Merge branch 'doself2-vectorisation' into dopair2-vectorisation

Conflicts:
	src/hydro/Gadget2/hydro_iact.h
	src/runner_doiact_vec.c
	src/runner_doiact_vec.h
	tests/test125cells.c
parents f4f24a0a 8193736e
......@@ -34,12 +34,17 @@ examples/*/*/*.txt
examples/*/*/used_parameters.yml
examples/*/gravity_checks_*.dat
tests/testPair
tests/testActivePair
tests/brute_force_periodic_BC_standard.dat
tests/swift_periodic_BC_standard.dat
tests/brute_force_periodic_BC_pertrubed.dat
tests/swift_periodic_BC_perturbed.dat
tests/brute_force_standard.dat
tests/swift_dopair_standard.dat
tests/brute_force_perturbed.dat
tests/swift_dopair_perturbed.dat
tests/test27cells
tests/testPeriodicBC
tests/test125cells
tests/brute_force_27_standard.dat
tests/swift_dopair_27_standard.dat
......@@ -49,6 +54,16 @@ tests/brute_force_125_standard.dat
tests/swift_dopair_125_standard.dat
tests/brute_force_125_perturbed.dat
tests/swift_dopair_125_perturbed.dat
tests/brute_force_active.dat
tests/brute_force_periodic_BC_perturbed.dat
tests/swift_dopair_active.dat
tests/test_nonsym_density_serial.dat
tests/test_nonsym_density_vec.dat
tests/test_nonsym_force_serial.dat
tests/test_nonsym_density_1_vec.dat
tests/test_nonsym_density_2_vec.dat
tests/test_nonsym_force_1_vec.dat
tests/test_nonsym_force_2_vec.dat
tests/testGreetings
tests/testReading
tests/input.hdf5
......@@ -64,12 +79,12 @@ tests/testMaths
tests/testThreadpool
tests/testParser
tests/parser_output.yml
tests/testPeriodicBC.sh
tests/testPeriodicBCPerturbed.sh
tests/test27cells.sh
tests/test27cellsPerturbed.sh
tests/test125cells.sh
tests/test125cellsPerturbed.sh
tests/testPair.sh
tests/testPairPerturbed.sh
tests/testParser.sh
tests/testReading.sh
tests/testAdiabaticIndex
......@@ -95,6 +110,9 @@ theory/paper_pasc/pasc_paper.pdf
theory/Multipoles/fmm.pdf
theory/Multipoles/fmm_standalone.pdf
theory/Multipoles/potential.pdf
theory/Multipoles/potential_long.pdf
theory/Multipoles/potential_short.pdf
theory/Multipoles/force_short.pdf
m4/libtool.m4
m4/ltoptions.m4
......
......@@ -15,28 +15,31 @@ Usage: swift [OPTION]... PARAMFILE
swift_mpi [OPTION]... PARAMFILE
Valid options are:
-a Pin runners using processor affinity.
-c Run with cosmological time integration.
-C Run with cooling.
-d Dry run. Read the parameter file, allocate memory but does not read
the particles from ICs and exit before the start of time integration.
Allows user to check validy of parameter and IC files as well as memory limits.
-D Always drift all particles even the ones far from active particles. This emulates
Gadget-[23] and GIZMO's default behaviours.
-e Enable floating-point exceptions (debugging mode).
-f {int} Overwrite the CPU frequency (Hz) to be used for time measurements.
-g Run with an external gravitational potential.
-G Run with self-gravity.
-M Reconstruct the multipoles every time-step.
-n {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop.
-s Run with hydrodynamics.
-S Run with stars.
-t {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
-T Print timers every time-step.
-v [12] Increase the level of verbosity
1: MPI-rank 0 writes
2: All MPI-ranks write
-y {int} Time-step frequency at which task graphs are dumped.
-h Print this help message and exit.
-a Pin runners using processor affinity.
-c Run with cosmological time integration.
-C Run with cooling.
-d Dry run. Read the parameter file, allocate memory but does not read
the particles from ICs and exit before the start of time integration.
Allows user to check validy of parameter and IC files as well as memory limits.
-D Always drift all particles even the ones far from active particles. This emulates
Gadget-[23] and GIZMO's default behaviours.
-e Enable floating-point exceptions (debugging mode).
-f {int} Overwrite the CPU frequency (Hz) to be used for time measurements.
-g Run with an external gravitational potential.
-G Run with self-gravity.
-M Reconstruct the multipoles every time-step.
-n {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop.
-P {sec:par:val} Set parameter value and overwrites values read from the parameters file. Can be used more than once.
-s Run with hydrodynamics.
-S Run with stars.
-t {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
-T Print timers every time-step.
-v [12] Increase the level of verbosity:
1: MPI-rank 0 writes,
2: All MPI-ranks write.
-y {int} Time-step frequency at which task graphs are dumped.
-Y {int} Time-step frequency at which threadpool tasks are dumped.
-h Print this help message and exit.
See the file parameter_example.yml for an example of parameter file.
See the file examples/parameter_example.yml for an example of parameter file.
......@@ -16,7 +16,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Init the project.
AC_INIT([SWIFT],[0.5.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
AC_INIT([SWIFT],[0.6.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
swift_config_flags="$*"
# Need to define this, instead of using fifth argument of AC_INIT, until 2.64.
......@@ -189,6 +189,19 @@ if test "$enable_task_debugging" = "yes"; then
AC_DEFINE([SWIFT_DEBUG_TASKS],1,[Enable task debugging])
fi
# Check if threadpool debugging is on.
AC_ARG_ENABLE([threadpool-debugging],
[AS_HELP_STRING([--enable-threadpool-debugging],
[Store threadpool mapper timing information and generate threadpool dump files @<:@yes/no@:>@]
)],
[enable_threadpool_debugging="$enableval"],
[enable_threadpool_debugging="no"]
)
if test "$enable_threadpool_debugging" = "yes"; then
AC_DEFINE([SWIFT_DEBUG_THREADPOOL],1,[Enable threadpool debugging])
LDFLAGS="$LDFLAGS -rdynamic"
fi
# Check if the general timers are switched on.
AC_ARG_ENABLE([timers],
[AS_HELP_STRING([--enable-timers],
......@@ -829,10 +842,10 @@ esac
# Gravity multipole order
AC_ARG_WITH([multipole-order],
[AS_HELP_STRING([--with-multipole-order=<order>],
[order of the multipole and gravitational field expansion @<:@ default: 3@:>@]
[order of the multipole and gravitational field expansion @<:@ default: 4@:>@]
)],
[with_multipole_order="$withval"],
[with_multipole_order="3"]
[with_multipole_order="4"]
)
AC_DEFINE_UNQUOTED([SELF_GRAVITY_MULTIPOLE_ORDER], [$with_multipole_order], [Multipole order])
......@@ -848,19 +861,31 @@ AM_CONDITIONAL([HAVE_DOXYGEN], [test "$ac_cv_path_ac_pt_DX_DOXYGEN" != ""])
# Handle .in files.
AC_CONFIG_FILES([Makefile src/Makefile examples/Makefile doc/Makefile doc/Doxyfile tests/Makefile])
AC_CONFIG_FILES([tests/testReading.sh], [chmod +x tests/testReading.sh])
AC_CONFIG_FILES([tests/testPair.sh], [chmod +x tests/testPair.sh])
AC_CONFIG_FILES([tests/testPairPerturbed.sh], [chmod +x tests/testPairPerturbed.sh])
AC_CONFIG_FILES([tests/testActivePair.sh], [chmod +x tests/testActivePair.sh])
AC_CONFIG_FILES([tests/test27cells.sh], [chmod +x tests/test27cells.sh])
AC_CONFIG_FILES([tests/test27cellsPerturbed.sh], [chmod +x tests/test27cellsPerturbed.sh])
AC_CONFIG_FILES([tests/test125cells.sh], [chmod +x tests/test125cells.sh])
AC_CONFIG_FILES([tests/test125cellsPerturbed.sh], [chmod +x tests/test125cellsPerturbed.sh])
AC_CONFIG_FILES([tests/testPeriodicBC.sh], [chmod +x tests/testPeriodicBC.sh])
AC_CONFIG_FILES([tests/testPeriodicBCPerturbed.sh], [chmod +x tests/testPeriodicBCPerturbed.sh])
AC_CONFIG_FILES([tests/testInteractions.sh], [chmod +x tests/testInteractions.sh])
AC_CONFIG_FILES([tests/testParser.sh], [chmod +x tests/testParser.sh])
# Save the compilation options
AC_DEFINE_UNQUOTED([SWIFT_CONFIG_FLAGS],["$swift_config_flags"],[Flags passed to configure])
# Make sure the latest git revision string gets included
touch src/version.c
# Generate output.
AC_OUTPUT
# Report general configuration.
AC_MSG_RESULT([
AC_MSG_RESULT([
------- Summary --------
$PACKAGE_NAME v.$PACKAGE_VERSION
Compiler : $CC
- vendor : $ax_cv_c_compiler_vendor
- version : $ax_cv_c_compiler_version
......@@ -887,14 +912,10 @@ AC_MSG_RESULT([
Multipole order : $with_multipole_order
No gravity below ID : $no_gravity_below_id
Individual timers : $enable_timers
Task debugging : $enable_task_debugging
Debugging checks : $enable_debugging_checks
Gravity checks : $gravity_force_checks
])
# Make sure the latest git revision string gets included
touch src/version.c
Individual timers : $enable_timers
Task debugging : $enable_task_debugging
Threadpool debugging : $enable_threadpool_debugging
Debugging checks : $enable_debugging_checks
Gravity checks : $gravity_force_checks
# Generate output.
AC_OUTPUT
------------------------])
......@@ -1988,6 +1988,9 @@ INCLUDE_FILE_PATTERNS =
# This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
PREDEFINED = "__attribute__(x)= "
PREDEFINED += HAVE_HDF5
PREDEFINED += WITH_MPI
PREDEFINED += WITH_VECTORIZATION
# If the MACRO_EXPANSION and EXPAND_ONLY_PREDEF tags are set to YES then this
# tag can be used to specify a list of macro names that should be expanded. The
......
......@@ -34,7 +34,7 @@ import sys
stats_filename = "./energy.txt"
# First snapshot
snap_filename = "coolingBox_000.hdf5"
snap_filename = "coolingBox_0000.hdf5"
# Some constants in cgs units
k_b = 1.38E-16 #boltzmann
......@@ -104,7 +104,7 @@ print "Cooling time:", cooling_time_cgs, "[s]"
u_snapshots_cgs = zeros(25)
t_snapshots_cgs = zeros(25)
for i in range(25):
snap = h5.File("coolingBox_%0.3d.hdf5"%i,'r')
snap = h5.File("coolingBox_%0.4d.hdf5"%i,'r')
u_snapshots_cgs[i] = sum(snap["/PartType0/InternalEnergy"][:] * snap["/PartType0/Masses"][:]) / total_mass[0] * unit_length**2 / (unit_time)**2
t_snapshots_cgs[i] = snap["/Header"].attrs["Time"] * unit_time
......
......@@ -20,7 +20,7 @@ H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
#read some header/parameter information from the first snapshot
filename = "Hydrostatic_000.hdf5"
filename = "Hydrostatic_0000.hdf5"
f = h5.File(filename,'r')
params = f["Parameters"]
unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
......@@ -39,7 +39,7 @@ M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
for i in range(n_snaps):
filename = "Hydrostatic_%03d.hdf5" %i
filename = "Hydrostatic_%04d.hdf5" %i
f = h5.File(filename,'r')
coords_dset = f["PartType0/Coordinates"]
coords = np.array(coords_dset)
......
......@@ -38,7 +38,7 @@ H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
#read some header/parameter information from the first snapshot
filename = "Hydrostatic_000.hdf5"
filename = "Hydrostatic_0000.hdf5"
f = h5.File(filename,'r')
params = f["Parameters"]
unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
......@@ -57,7 +57,7 @@ M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
for i in range(n_snaps):
filename = "Hydrostatic_%03d.hdf5" %i
filename = "Hydrostatic_%04d.hdf5" %i
f = h5.File(filename,'r')
coords_dset = f["PartType0/Coordinates"]
coords = np.array(coords_dset)
......
......@@ -17,7 +17,7 @@ H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
#read some header/parameter information from the first snapshot
filename = "CoolingHalo_000.hdf5"
filename = "CoolingHalo_0000.hdf5"
f = h5.File(filename,'r')
params = f["Parameters"]
unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
......@@ -41,7 +41,7 @@ time_array_cgs = []
for i in range(n_snaps):
filename = "CoolingHalo_%03d.hdf5" %i
filename = "CoolingHalo_%04d.hdf5" %i
f = h5.File(filename,'r')
coords_dset = f["PartType0/Coordinates"]
coords = np.array(coords_dset)
......
......@@ -39,7 +39,7 @@ H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
#read some header/parameter information from the first snapshot
filename = "CoolingHalo_000.hdf5"
filename = "CoolingHalo_0000.hdf5"
f = h5.File(filename,'r')
params = f["Parameters"]
unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
......@@ -58,7 +58,7 @@ M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
for i in range(n_snaps):
filename = "CoolingHalo_%03d.hdf5" %i
filename = "CoolingHalo_%04d.hdf5" %i
f = h5.File(filename,'r')
coords_dset = f["PartType0/Coordinates"]
coords = np.array(coords_dset)
......
......@@ -21,7 +21,7 @@ H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
#read some header/parameter information from the first snapshot
filename = "CoolingHalo_000.hdf5"
filename = "CoolingHalo_0000.hdf5"
f = h5.File(filename,'r')
params = f["Parameters"]
unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
......@@ -42,7 +42,7 @@ M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
for i in range(n_snaps):
filename = "CoolingHalo_%03d.hdf5" %i
filename = "CoolingHalo_%04d.hdf5" %i
f = h5.File(filename,'r')
coords_dset = f["PartType0/Coordinates"]
coords = np.array(coords_dset)
......
......@@ -39,7 +39,7 @@ H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
#read some header/parameter information from the first snapshot
filename = "CoolingHalo_000.hdf5"
filename = "CoolingHalo_0000.hdf5"
f = h5.File(filename,'r')
params = f["Parameters"]
unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
......@@ -60,7 +60,7 @@ M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
for i in range(n_snaps):
filename = "CoolingHalo_%03d.hdf5" %i
filename = "CoolingHalo_%04d.hdf5" %i
f = h5.File(filename,'r')
coords_dset = f["PartType0/Coordinates"]
coords = np.array(coords_dset)
......
......@@ -20,7 +20,7 @@ H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
#read some header/parameter information from the first snapshot
filename = "CoolingHalo_000.hdf5"
filename = "CoolingHalo_0000.hdf5"
f = h5.File(filename,'r')
params = f["Parameters"]
unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
......@@ -44,7 +44,7 @@ time_array_cgs = []
for i in range(n_snaps):
filename = "CoolingHalo_%03d.hdf5" %i
filename = "CoolingHalo_%04d.hdf5" %i
f = h5.File(filename,'r')
coords_dset = f["PartType0/Coordinates"]
coords = np.array(coords_dset)
......
......@@ -39,7 +39,7 @@ H_0_cgs = 100. * h * KM_PER_SEC_IN_CGS / (1.0e6 * PARSEC_IN_CGS)
#read some header/parameter information from the first snapshot
filename = "CoolingHalo_000.hdf5"
filename = "CoolingHalo_0000.hdf5"
f = h5.File(filename,'r')
params = f["Parameters"]
unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
......@@ -58,7 +58,7 @@ M_vir_cgs = r_vir_cgs * v_c_cgs**2 / CONST_G_CGS
for i in range(n_snaps):
filename = "CoolingHalo_%03d.hdf5" %i
filename = "CoolingHalo_%04d.hdf5" %i
f = h5.File(filename,'r')
coords_dset = f["PartType0/Coordinates"]
coords = np.array(coords_dset)
......
......@@ -18,3 +18,5 @@ output to 'Disc-Patch-dynamic.hdf5'. These are now the ICs for the actual test.
When running SWIFT with the parameters from 'disc-patch.yml' and an
ideal gas EoS on these ICs the disc should stay in equilibrium.
The solution can be checked using the 'plotSolution.py' script.
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.9885e33 # Grams
UnitLength_in_cgs: 3.0856776e18 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitMass_in_cgs: 1.9885e33 # Grams
UnitLength_in_cgs: 3.08567758149e18 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
......@@ -11,17 +11,17 @@ TimeIntegration:
time_begin: 0 # The starting time of the simulation (in internal units).
time_end: 968. # The end time of the simulation (in internal units).
dt_min: 1e-4 # The minimal time-step size of the simulation (in internal units).
dt_max: 1. # The maximal time-step size of the simulation (in internal units).
dt_max: 10. # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1 # Time between statistics output
delta_time: 12. # Time between statistics output
# Parameters governing the snapshots
Snapshots:
basename: Disc-Patch # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 12. # Time difference between consecutive outputs (in internal units)
basename: Disc-Patch # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 48. # Time difference between outputs (in internal units)
# Parameters for the hydrodynamics scheme
SPH:
......@@ -29,7 +29,7 @@ SPH:
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
max_smoothing_length: 70. # Maximal smoothing length allowed (in internal units).
h_max: 60. # Maximal smoothing length allowed (in internal units).
# Parameters related to the initial conditions
InitialConditions:
......@@ -39,6 +39,8 @@ InitialConditions:
DiscPatchPotential:
surface_density: 10.
scale_height: 100.
z_disc: 200.
x_disc: 400.
x_trunc: 300.
x_max: 350.
timestep_mult: 0.03
growth_time: 5.
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.9885e33 # Grams
UnitLength_in_cgs: 3.0856776e18 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitMass_in_cgs: 1.9885e33 # Grams
UnitLength_in_cgs: 3.08567758149e18 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
......@@ -11,17 +11,17 @@ TimeIntegration:
time_begin: 968 # The starting time of the simulation (in internal units).
time_end: 12000. # The end time of the simulation (in internal units).
dt_min: 1e-4 # The minimal time-step size of the simulation (in internal units).
dt_max: 1. # The maximal time-step size of the simulation (in internal units).
dt_max: 10. # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1 # Time between statistics output
delta_time: 24 # Time between statistics output
# Parameters governing the snapshots
Snapshots:
basename: Disc-Patch-dynamic # Common part of the name of output files
time_first: 968. # Time of the first output (in internal units)
delta_time: 24. # Time difference between consecutive outputs (in internal units)
basename: Disc-Patch-dynamic # Common part of the name of output files
time_first: 968. # Time of the first output (in internal units)
delta_time: 96. # Time difference between outputs (in internal units)
# Parameters for the hydrodynamics scheme
SPH:
......@@ -29,7 +29,7 @@ SPH:
delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
max_ghost_iterations: 30 # Maximal number of iterations allowed to converge towards the smoothing length.
max_smoothing_length: 70. # Maximal smoothing length allowed (in internal units).
h_max: 60. # Maximal smoothing length allowed (in internal units).
# Parameters related to the initial conditions
InitialConditions:
......@@ -39,5 +39,7 @@ InitialConditions:
DiscPatchPotential:
surface_density: 10.
scale_height: 100.
z_disc: 200.
x_disc: 400.
x_trunc: 300.
x_max: 380.
timestep_mult: 0.03
;
; test energy / angular momentum conservation of test problem
;
iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare final and initial energy
; set physical constants
@physunits
indir = './'
;basefile = 'Disc-Patch-dynamic_'
basefile = 'Disc-Patch_'
; set properties of potential
uL = phys.pc ; unit of length
uM = phys.msun ; unit of mass
uV = 1d5 ; unit of velocity
; properties of patch
surface_density = 100. ; surface density of all mass, which generates the gravitational potential
scale_height = 100.
z_disk = 200. ;
fgas = 0.1 ; gas fraction
gamma = 5./3.
; derived units
constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
pcentre = [0.,0.,z_disk] * pc / uL
utherm = !pi * constG * surface_density * scale_height / (gamma-1.)
temp = (utherm*uV^2)*phys.m_h/phys.kb
soundspeed = sqrt(gamma * (gamma-1.) * utherm)
t_dyn = sqrt(scale_height / (constG * surface_density))
rho0 = fgas*(surface_density)/(2.*scale_height)
print,' dynamical time = ',t_dyn,' = ',t_dyn*UL/uV/(1d6*phys.yr),' Myr'
print,' thermal energy per unit mass = ',utherm
print,' central density = ',rho0,' = ',rho0*uM/uL^3/m_h,' particles/cm^3'
print,' central temperature = ',temp
lambda = 2 * !pi * phys.G^1.5 * (scale_height*uL)^1.5 * (surface_density * uM/uL^2)^0.5 * phys.m_h^2 / (gamma-1) / fgas
print,' lambda = ',lambda
stop
;
infile = indir + basefile + '*'
spawn,'ls -1 '+infile,res
nfiles = n_elements(res)
; choose: calculate change of energy and Lz, comparing first and last
; snapshots for all particles, or do so for a subset
; compare all
ifile = 0
inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
id = h5rd(inf,'PartType0/ParticleIDs')
nfollow = n_elements(id)
; compute anlytic profile
nbins = 100
zbins = findgen(nbins)/float(nbins-1) * 2 * scale_height
rbins = (surface_density/(2.*scale_height)) / cosh(abs(zbins)/scale_height)^2
; plot analytic profile
wset,0
plot,[0],[0],xr=[0,2*scale_height],yr=[0,max(rbins)],/nodata,xtitle='|z|',ytitle=textoidl('\rho')
oplot,zbins,rbins,color=blue
ifile = 0
nskip = nfiles - 1
isave = 0
nplot = 8192 ; randomly plot particles
color = floor(findgen(nfiles)/float(nfiles-1)*255)
;for ifile=0,nfiles-1,nskip do begin
tsave = [0.]
toplot = [1,nfiles-1]
for iplot=0,1 do begin
ifile = toplot[iplot]
inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
time = h5ra(inf, 'Header','Time')
tsave = [tsave, time]
print,' time= ',time
p = h5rd(inf,'PartType0/Coordinates')
v = h5rd(inf,'PartType0/Velocities')
id = h5rd(inf,'PartType0/ParticleIDs')
rho = h5rd(inf,'PartType0/Density')
h = h5rd(inf,'PartType0/SmoothingLength')
utherm = h5rd(inf,'PartType0/InternalEnergy')
indx = sort(id)
; substract disk centre
for ic=0,2 do p[ic,*]=p[ic,*] - pcentre[ic]
;; ; if you want to sort particles by ID
;; id = id[indx]
;; rho = rho[indx]
;; utherm = utherm[indx]
;; h = h[indx]
;; for ic=0,2 do begin
;; tmp = reform(p[ic,*]) & p[ic,*] = tmp[indx]
;; tmp = reform(v[ic,*]) & v[ic,*] = tmp[indx]
;; endfor
ip = floor(randomu(ifile+1,nplot)*n_elements(rho))
color = red
if(ifile eq 1) then begin
color=black
endif else begin
color=red
endelse
oplot,abs(p[2,ip]), rho[ip], psym=3, color=color
isave = isave + 1
endfor
; time in units of dynamical time