Merge branch 'master' into FOF

Conflicts:
	examples/main.c
	src/Makefile.am
	src/engine.c
	src/engine.h

.clang-format

@@ -2,4 +2,5 @@
 Language:        Cpp
 BasedOnStyle:  Google
 KeepEmptyLinesAtTheStartOfBlocks: true
+PenaltyBreakAssignment: 2
 ...

.gitignore

@@ -42,6 +42,11 @@ examples/*/*/*.hdf5
 examples/*/snapshots*
 examples/*/restart/*
 examples/*/*/used_parameters.yml
+examples/*/err_file*
+examples/*/out_file*
+examples/*/stf_output*
+examples/*/stf_ouput*
+examples/*/log*
 examples/*/*/unused_parameters.yml
 examples/*/*.mpg
 examples/*/gravity_checks_*.dat
@@ -106,6 +111,8 @@ tests/test125cellsPerturbed.sh
 tests/testParser.sh
 tests/testReading.sh
 tests/testSelectOutput.sh
+tests/unused_parser_output.yml
+tests/used_parser_output.yml
 tests/testAdiabaticIndex
 tests/testRiemannExact
 tests/testRiemannTRRS
@@ -124,6 +131,7 @@ tests/testEOS
 tests/testEOS*.txt
 tests/testEOS*.png
 tests/testUtilities
+tests/testCbrt
 
 theory/latex/swift.pdf
 theory/SPH/Kernels/kernels.pdf
@@ -157,6 +165,8 @@ m4/lt~obsolete.m4
 /stamp-h1
 /test-driver
 
+src/equation_of_state/planetary/*.txt
+
 # Intel compiler optimization reports
 *.optrpt
 
@@ -308,3 +318,6 @@ sympy-plots-for-*.tex/
 
 #ctags
 *tags
+
+# vim
+*.swp

AUTHORS

@@ -2,12 +2,13 @@ Pedro Gonnet            gonnet@google.com
 Matthieu Schaller       matthieu.schaller@durham.ac.uk
 Aidan Chalk             aidan.chalk@durham.ac.uk
 Peter W. Draper         p.w.draper@durham.ac.uk
-Bert Vandenbrouck       bert.vandenbroucke@gmail.com
+Bert Vandenbroucke      bert.vandenbroucke@gmail.com
 James S. Willis 	james.s.willis@durham.ac.uk
 John A. Regan 		john.a.regan@durham.ac.uk
 Angus Lepper		angus.lepper@ed.ac.uk
 Tom Theuns 		tom.theuns@durham.ac.uk
 Richard G. Bower	r.g.bower@durham.ac.uk
 Stefan Arridge		stefan.arridge@durham.ac.uk
-Massimiliano Culpo	massimiliano.culpo@googlemail.com
+Josh Borrow             joshua.borrow@durham.ac.uk
+Loic Hausammann		loic.hausammann@epfl.ch
 Yves Revaz   		yves.revaz@epfl.ch

CONTRIBUTING.md

 The SWIFT source code is using a variation of the 'Google' formatting style. 
-The script 'format.sh' in the root directory applies the clang-format-3.8
+The script 'format.sh' in the root directory applies the clang-format-5.0
 tool with our style choices to all the SWIFT C source file. Please apply 
 the formatting script to the files before submitting a merge request.
 
@@ -9,4 +9,4 @@ check` in the root directory. Please check that the test suite still
 runs with your changes applied before submitting a merge request and add 
 relevant unit tests probing the correctness of new modules. An example of how
 to add a test to the suite can be found by considering the tests/testGreeting 
-case.
+case.
\ No newline at end of file

INSTALL.swift

@@ -182,6 +182,6 @@ before you can build it.
                              ==================
 
 The SWIFT source code uses a variation of 'Google' style. The script
-'format.sh' in the root directory applies the clang-format-3.8 tool with our
+'format.sh' in the root directory applies the clang-format-5.0 tool with our
 style choices to all the SWIFT C source file. Please apply the formatting
 script to the files before submitting a merge request.

README

@@ -29,7 +29,9 @@ Valid options are:
   -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.
+  -o          {str} Generate a default output parameter file.
   -P  {sec:par:val} Set parameter value and overwrites values read from the parameters file. Can be used more than once.
+  -r                Continue using restart files.
   -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.
@@ -37,6 +39,7 @@ Valid options are:
   -v           [12] Increase the level of verbosity:
                     1: MPI-rank 0 writes,
                     2: All MPI-ranks write.
+  -x                Run with structure finding.
   -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.

README.md

@@ -32,7 +32,7 @@ Contribution Guidelines
 -----------------------
 
 The SWIFT source code uses a variation of the 'Google' formatting style.
-The script 'format.sh' in the root directory applies the clang-format-3.8
+The script 'format.sh' in the root directory applies the clang-format-5.0
 tool with our style choices to all the SWIFT C source file. Please apply
 the formatting script to the files before submitting a pull request.
 
@@ -77,7 +77,9 @@ Valid options are:
   -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.
+  -o          {str} Generate a default output parameter file.
   -P  {sec:par:val} Set parameter value and overwrites values read from the parameters file. Can be used more than once.
+  -r                Continue using restart files.
   -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.
@@ -85,6 +87,7 @@ Valid options are:
   -v           [12] Increase the level of verbosity:
                     1: MPI-rank 0 writes,
                     2: All MPI-ranks write.
+  -x                Run with structure finding.
   -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.

configure.ac

@@ -38,47 +38,6 @@ AX_CHECK_ENABLE_DEBUG
 AC_PROG_CC
 AM_PROG_CC_C_O
 
-
-# Master subgrid options
-# If you add a restriction (e.g. no cooling, chemistry or hydro)
-# you will need to check for overwrite after reading the additional options.
-# As an example for this, see the call to AC_ARG_WITH for cooling.
-AC_ARG_WITH([subgrid],
-	[AS_HELP_STRING([--with-subgrid=<subgrid>],
-		[Master switch for subgrid methods. Inexperienced user should start from here @<:@none, GEAR, EAGLE default: none@:>@]
-	)],
-	[with_subgrid="$withval"],
-	[with_subgrid=none]
-)
-
-# Default values
-with_subgrid_cooling=none
-with_subgrid_chemistry=none
-with_subgrid_hydro=none
-
-case "$with_subgrid" in
-   yes)
-      AC_MSG_ERROR([Invalid option. A subgrid model must be chosen.])
-   ;;
-   none)
-   ;;
-   GEAR)
-	with_subgrid_cooling=grackle
-	with_subgrid_chemistry=GEAR
-	with_subgrid_hydro=gadget2
-   ;;
-   EAGLE)
-	with_subgrid_cooling=EAGLE
-	with_subgrid_chemistry=EAGLE
-	with_subgrid_hydro=gadget2
-   ;;
-   *)
-      AC_MSG_ERROR([Unknown subgrid choice: $with_subgrid])
-   ;;
-esac
-
-
-
 # If debug is selected then we also define SWIFT_DEVELOP_MODE to control
 # any developer code options.
 if test "x$ax_enable_debug" != "xno"; then
@@ -273,6 +232,18 @@ if test "$enable_debugging_checks" = "yes"; then
    AC_DEFINE([SWIFT_DEBUG_CHECKS],1,[Enable expensive debugging])
 fi
 
+# Check if using our custom icbrtf is enalbled.
+AC_ARG_ENABLE([custom-icbrtf],
+   [AS_HELP_STRING([--enable-custom-icbrtf],
+     [Use SWIFT's custom icbrtf function instead of the system cbrtf @<:@yes/no@:>@]
+   )],
+   [enable_custom_icbrtf="$enableval"],
+   [enable_custom_icbrtf="no"]
+)
+if test "$enable_custom_icbrtf" = "yes"; then
+   AC_DEFINE([WITH_ICBRTF],1,[Enable custom icbrtf])
+fi
+
 # Check whether we want to default to naive cell interactions
 AC_ARG_ENABLE([naive-interactions],
    [AS_HELP_STRING([--enable-naive-interactions],
@@ -549,42 +520,6 @@ AH_VERBATIM([__STDC_FORMAT_MACROS],
 #define __STDC_FORMAT_MACROS 1
 #endif])
 
-# Check for grackle.
-have_grackle="no"
-AC_ARG_WITH([grackle],
-    [AS_HELP_STRING([--with-grackle=PATH],
-       [root directory where grackle is installed @<:@yes/no@:>@]
-    )],
-    [with_grackle="$withval"],
-    [with_grackle="no"]
-)
-if test "x$with_grackle" != "xno"; then
-   AC_PROG_FC
-   AC_FC_LIBRARY_LDFLAGS
-   if test "x$with_grackle" != "xyes" -a "x$with_grackle" != "x"; then
-      GRACKLE_LIBS="-L$with_grackle/lib -lgrackle"
-      GRACKLE_INCS="-I$with_grackle/include"
-   else
-      GRACKLE_LIBS="-lgrackle"
-      GRACKLE_INCS=""
-   fi
-
-   have_grackle="yes"
-
-   AC_CHECK_LIB(
-      [grackle],
-      [initialize_chemistry_data],
-      [AC_DEFINE([HAVE_GRACKLE],1,[The GRACKLE library appears to be present.])
-        AC_DEFINE([CONFIG_BFLOAT_8],1,[Use doubles in grackle])
-      ],
-      [AC_MSG_ERROR(Cannot find grackle library!)],
-      [$GRACKLE_LIBS $GRACKLE_INCS $FCLIBS]
-   )
-fi
-AC_SUBST([GRACKLE_LIBS])
-AC_SUBST([GRACKLE_INCS])
-AM_CONDITIONAL([HAVEGRACKLE],[test -n "$GRACKLE_LIBS"])
-
 # Check for FFTW. We test for this in the standard directories by default,
 # and only disable if using --with-fftw=no or --without-fftw. When a value
 # is given GSL must be found.
@@ -641,7 +576,7 @@ if test "x$with_profiler" != "xno"; then
       proflibs="-lprofiler"
    fi
    AC_CHECK_LIB([profiler],[ProfilerFlush],
-    [have_profiler="yes" 
+    [have_profiler="yes"
       AC_DEFINE([WITH_PROFILER],1,[Link against the gperftools profiling library.])],
     [have_profiler="no"], $proflibs)
 
@@ -836,6 +771,38 @@ if test "$with_hdf5" = "yes"; then
 fi
 AM_CONDITIONAL([HAVEPARALLELHDF5],[test "$have_parallel_hdf5" = "yes"])
 
+# Check for VELOCIraptor.
+have_velociraptor="no"
+AC_ARG_WITH([velociraptor],
+    [AS_HELP_STRING([--with-velociraptor=PATH],
+       [Directory where velociraptor library exists @<:@yes/no@:>@]
+    )],
+    [with_velociraptor="$withval"],
+    [with_velociraptor="no"]
+)
+if test "x$with_velociraptor" != "xno"; then
+   AC_PROG_FC
+   AC_FC_LIBRARY_LDFLAGS
+   if test "x$with_velociraptor" != "xyes" -a "x$with_velociraptor" != "x"; then
+      VELOCIRAPTOR_LIBS="-L$with_velociraptor -lstf -lstdc++ -lhdf5_cpp"
+      CFLAGS="$CFLAGS -fopenmp"
+   else
+      VELOCIRAPTOR_LIBS=""
+   fi
+
+   have_velociraptor="yes"
+
+   AC_CHECK_LIB(
+      [stf],
+      [InitVelociraptor],
+      [AC_DEFINE([HAVE_VELOCIRAPTOR],1,[The VELOCIraptor library appears to be present.])],
+      [AC_MSG_ERROR(Cannot find VELOCIraptor library at $with_velociraptor)],
+      [$VELOCIRAPTOR_LIBS $HDF5_LDFLAGS $HDF5_LIBS $GSL_LIBS]
+   )
+fi
+AC_SUBST([VELOCIRAPTOR_LIBS])
+AM_CONDITIONAL([HAVEVELOCIRAPTOR],[test -n "$VELOCIRAPTOR_LIBS"])
+
 # Check for floating-point execeptions
 AC_CHECK_FUNC(feenableexcept, AC_DEFINE([HAVE_FE_ENABLE_EXCEPT],[1],
     [Defined if the floating-point exception can be enabled using non-standard GNU functions.]))
@@ -944,10 +911,69 @@ fi
 
 # Various package configuration options.
 
+# Master subgrid options
+# If you add a restriction (e.g. no cooling, chemistry or hydro)
+# you will need to check for overwrite after reading the additional options.
+# As an example for this, see the call to AC_ARG_WITH for cooling.
+AC_ARG_WITH([subgrid],
+	[AS_HELP_STRING([--with-subgrid=<subgrid>],
+		[Master switch for subgrid methods. Inexperienced user should start from here @<:@none, GEAR, EAGLE default: none@:>@]
+	)],
+	[with_subgrid="$withval"],
+	[with_subgrid=none]
+)
+
+# Default values
+with_subgrid_cooling=none
+with_subgrid_chemistry=none
+with_subgrid_hydro=none
+
+case "$with_subgrid" in
+   yes)
+      AC_MSG_ERROR([Invalid option. A subgrid model must be chosen.])
+   ;;
+   none)
+   ;;
+   GEAR)
+	with_subgrid_cooling=grackle
+	with_subgrid_chemistry=GEAR
+	with_subgrid_hydro=gadget2
+   ;;
+   EAGLE)
+	with_subgrid_cooling=EAGLE
+	with_subgrid_chemistry=EAGLE
+	with_subgrid_hydro=gadget2
+   ;;
+   *)
+      AC_MSG_ERROR([Unknown subgrid choice: $with_subgrid])
+   ;;
+esac
+
+# Gravity scheme.
+AC_ARG_WITH([gravity],
+   [AS_HELP_STRING([--with-gravity=<scheme>],
+      [Gravity scheme to use @<:@default, with-potential, default: default@:>@]
+   )],
+   [with_gravity="$withval"],
+   [with_gravity="default"]
+)
+
+case "$with_gravity" in
+   with-potential)
+      AC_DEFINE([POTENTIAL_GRAVITY], [1], [Gravity scheme with potential calculation])
+   ;;
+   default)
+      AC_DEFINE([DEFAULT_GRAVITY], [1], [Default gravity scheme])
+   ;;
+   *)
+      AC_MSG_ERROR([Unknown gravity scheme: $with_gravity])
+   ;;
+esac
+
 # Hydro scheme.
 AC_ARG_WITH([hydro],
    [AS_HELP_STRING([--with-hydro=<scheme>],
-      [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, default, gizmo-mfv, gizmo-mfm, shadowfax, minimal-multi-mat, debug default: gadget2@:>@]
+      [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, debug default: gadget2@:>@]
    )],
    [with_hydro="$withval"],
    [with_hydro="gadget2"]
@@ -986,10 +1012,11 @@ case "$with_hydro" in
    shadowfax)
       AC_DEFINE([SHADOWFAX_SPH], [1], [Shadowfax SPH])
    ;;
-   minimal-multi-mat)
-      AC_DEFINE([MINIMAL_MULTI_MAT_SPH], [1], [Minimal Multiple Material SPH])
+   planetary)
+      AC_DEFINE([PLANETARY_SPH], [1], [Planetary SPH])
    ;;
 
+
    *)
       AC_MSG_ERROR([Unknown hydrodynamics scheme: $with_hydro])
    ;;
@@ -1148,6 +1175,42 @@ case "$with_riemann" in
    ;;
 esac
 
+# Check for grackle.
+have_grackle="no"
+AC_ARG_WITH([grackle],
+    [AS_HELP_STRING([--with-grackle=PATH],
+       [root directory where grackle is installed @<:@yes/no@:>@]
+    )],
+    [with_grackle="$withval"],
+    [with_grackle="no"]
+)
+if test "x$with_grackle" != "xno"; then
+   AC_PROG_FC
+   AC_FC_LIBRARY_LDFLAGS
+   if test "x$with_grackle" != "xyes" -a "x$with_grackle" != "x"; then
+      GRACKLE_LIBS="-L$with_grackle/lib -lgrackle"
+      GRACKLE_INCS="-I$with_grackle/include"
+   else
+      GRACKLE_LIBS="-lgrackle"
+      GRACKLE_INCS=""
+   fi
+
+   have_grackle="yes"
+
+   AC_CHECK_LIB(
+      [grackle],
+      [initialize_chemistry_data],
+      [AC_DEFINE([HAVE_GRACKLE],1,[The GRACKLE library appears to be present.])
+        AC_DEFINE([CONFIG_BFLOAT_8],1,[Use doubles in grackle])
+      ],
+      [AC_MSG_ERROR(Cannot find grackle library!)],
+      [$GRACKLE_LIBS $GRACKLE_INCS $FCLIBS]
+   )
+fi
+AC_SUBST([GRACKLE_LIBS])
+AC_SUBST([GRACKLE_INCS])
+AM_CONDITIONAL([HAVEGRACKLE],[test -n "$GRACKLE_LIBS"])
+
 #  Cooling function
 AC_ARG_WITH([cooling],
    [AS_HELP_STRING([--with-cooling=<function>],
@@ -1269,10 +1332,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: 5@:>@]
+      [order of the multipole and gravitational field expansion @<:@ default: 4@:>@]
    )],
    [with_multipole_order="$withval"],
-   [with_multipole_order="5"]
+   [with_multipole_order="4"]
 )
 AC_DEFINE_UNQUOTED([SELF_GRAVITY_MULTIPOLE_ORDER], [$with_multipole_order], [Multipole order])
 
@@ -1318,21 +1381,22 @@ AC_MSG_RESULT([
 
    $PACKAGE_NAME v.$PACKAGE_VERSION
 
-   Compiler           : $CC
-    - vendor          : $ax_cv_c_compiler_vendor
-    - version         : $ax_cv_c_compiler_version
-    - flags           : $CFLAGS
-   MPI enabled        : $enable_mpi
-   HDF5 enabled       : $with_hdf5
-    - parallel        : $have_parallel_hdf5
-   Metis enabled      : $have_metis
-   FFTW3 enabled      : $have_fftw
-   GSL enabled        : $have_gsl
-   libNUMA enabled    : $have_numa
-   GRACKLE enabled    : $have_grackle
-   Special allocators : $have_special_allocator
-   CPU profiler       : $have_profiler
-   Pthread barriers   : $have_pthread_barrier
+   Compiler             : $CC
+    - vendor            : $ax_cv_c_compiler_vendor
+    - version           : $ax_cv_c_compiler_version
+    - flags             : $CFLAGS
+   MPI enabled          : $enable_mpi
+   HDF5 enabled         : $with_hdf5
+    - parallel          : $have_parallel_hdf5
+   Metis enabled        : $have_metis
+   FFTW3 enabled        : $have_fftw
+   GSL enabled          : $have_gsl
+   libNUMA enabled      : $have_numa
+   GRACKLE enabled      : $have_grackle
+   Special allocators   : $have_special_allocator
+   CPU profiler         : $have_profiler
+   Pthread barriers     : $have_pthread_barrier
+   VELOCIraptor enabled : $have_velociraptor
 
    Hydro scheme       : $with_hydro
    Dimensionality     : $with_dimension
@@ -1340,12 +1404,14 @@ AC_MSG_RESULT([
    Equation of state  : $with_eos
    Adiabatic index    : $with_gamma
    Riemann solver     : $with_riemann
-   Cooling function   : $with_cooling
-   Chemistry          : $with_chemistry
 
-   External potential  : $with_potential
+   Gravity scheme      : $with_gravity
    Multipole order     : $with_multipole_order
    No gravity below ID : $no_gravity_below_id
+   External potential  : $with_potential
+
+   Cooling function   : $with_cooling
+   Chemistry          : $with_chemistry
 
    Individual timers     : $enable_timers
    Task debugging        : $enable_task_debugging
@@ -1354,5 +1420,6 @@ AC_MSG_RESULT([
    Interaction debugging : $enable_debug_interactions
    Naive interactions    : $enable_naive_interactions
    Gravity checks        : $gravity_force_checks
+   Custom icbrtf         : $enable_custom_icbrtf
 
  ------------------------])

doc/Doxyfile.in

@@ -766,8 +766,8 @@ INPUT		       += @top_srcdir@/src/stars/Default
 INPUT		       += @top_srcdir@/src/riemann
 INPUT		       += @top_srcdir@/src/potential/point_mass
 INPUT		       += @top_srcdir@/src/equation_of_state/ideal_gas
-INPUT		       += @top_srcdir@/src/cooling/none
-INPUT		       += @top_srcdir@/src/chemistry/none
+INPUT		       += @top_srcdir@/src/cooling/EAGLE
+INPUT		       += @top_srcdir@/src/chemistry/EAGLE
 
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
@@ -818,7 +818,7 @@ EXCLUDE_SYMLINKS       = NO
 # Note that the wildcards are matched against the file with absolute path, so to
 # exclude all test directories for example use the pattern */test/*
 
-EXCLUDE_PATTERNS       =
+EXCLUDE_PATTERNS       = *.md
 
 # The EXCLUDE_SYMBOLS tag can be used to specify one or more symbol names
 # (namespaces, classes, functions, etc.) that should be excluded from the
@@ -1957,7 +1957,7 @@ MACRO_EXPANSION        = YES
 # The default value is: NO.
 # This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
 
-EXPAND_ONLY_PREDEF     = YES
+EXPAND_ONLY_PREDEF     = NO
 
 # If the SEARCH_INCLUDES tag is set to YES, the include files in the
 # INCLUDE_PATH will be searched if a #include is found.
@@ -1991,8 +1991,10 @@ INCLUDE_FILE_PATTERNS  =
 
 PREDEFINED             = "__attribute__(x)= "
 PREDEFINED             += HAVE_HDF5
+PREDEFINED             += HAVE_FFTW
 PREDEFINED             += WITH_MPI
 PREDEFINED             += WITH_VECTORIZATION
+PREDEFINED             += __GNUC__
 
 # 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

doc/RTD/source/GettingStarted/parameter_file.rst

 .. Parameter File
    Loic Hausammann, 1 june 2018
 
+.. _Parameter_File_label:
+
 Parameter File
 ==============
 
@@ -10,6 +12,35 @@ parameters.  Each section in this file corresponds to a different option in
 SWIFT and are not always required depending on the configuration options and
 the run time parameters.
 
+Output List
+~~~~~~~~~~~
+
+In the sections ``Snapshots`` and ``Statistics``, you can specify the options ``output_list_on`` and ``output_list``  which receive an int and a filename.
+The ``output_list_on`` enable or not the output list and ``output_list`` is the filename containing the output times.
+With the file header, you can choose between writing redshifts, scale factors or times.
+
+Example of file containing with times (in internal units):
+::
+   # Time
+   0.5
+   1.5
+   3.0
+   12.5
+
+Example of file with scale factors:
+::
+   # Scale Factor
+   0.1
+   0.2
+   0.3
+
+Example of file with redshift: 
+::
+   # Redshift
+   20
+   15
+   10
+   5
 
 Output Selection
 ~~~~~~~~~~~~~~~~

doc/RTD/source/InitialConditions/index.rst

@@ -4,45 +4,63 @@
 Initial Conditions
 ==================
 
-To run anything more than examples from our suite, you will need to be able to 
-produce your own initial conditions for SWIFT. We use the same initial conditions
-as the popular GADGET-2 code, which uses the HDF5 file format.
+To run anything more than examples from our suite, you will need to be able to
+produce your own initial conditions for SWIFT. We use the same initial
+conditions format as the popular `GADGET-2
+<https://wwwmpa.mpa-garching.mpg.de/~volker/gadget/>`_ code, which uses HDF5 for
+its type 3 format. Note that we do not support the GADGET-2 types 1 and 2
+formats.
+
+The original GADGET-2 file format only contains 2 types of particles: gas
+particles and 5 sorts of collisionless particles that allow users to run with 5
+separate particle masses and softenings. In SWIFT, we expand on this by using
+two of these types for stars and black holes.
+
+GADGET-2 can have initial conditions split over many files. This allow multiple
+ones to be read in parallel and is the only way the code can handle more than
+2^31 particles. This limitation is not in place in SWIFT. A single file can
+contain any number of particles (well... up to 2^64...) and the file is read in
+parallel by HDF5 when running on more than one compute node.
 
 As the original documentation for the GADGET-2 initial conditions format is
-quite sparse, we lay out here all of the necessary components. If you are generating
-your initial conditions from python, we recommend you use the h5py package. We
-provide a writing wrapper for this for our initial conditions in
+quite sparse, we lay out here all of the necessary components. If you are
+generating your initial conditions from python, we recommend you use the h5py
+package. We provide a writing wrapper for this for our initial conditions in
 ``examples/KeplerianRing/write_gadget.py``.
 
-You can find out more about the HDF5 format on their webpages, here:
-https://support.hdfgroup.org/HDF5/doc/H5.intro.html
+You can find out more about the HDF5 format on their `webpages
+<https://support.hdfgroup.org/HDF5/doc/H5.intro.html>`_.
 
 
 Structure of the File
 ---------------------
 
-There are several groups that contain 'auxilliary' information, such as ``Header``.
-Particle data is placed in groups that signify particle type.
-
-+---------------------+------------------------+
-| Group Name          | Physical Particle Type |
-+=====================+========================+
-| ``PartType0``       | Gas                    |
-+---------------------+------------------------+
-| ``PartType1``       | Dark Matter            |
-+---------------------+------------------------+
-| ``PartType2``       | Ignored                |
-+---------------------+------------------------+
-| ``PartType3``       | Ignored                |
-+---------------------+------------------------+
-| ``PartType4``       | Stars                  |
-+---------------------+------------------------+
-| ``PartType5``       | Black Holes            |
-+---------------------+------------------------+
-
-Currently, not all of these particle types are included in SWIFT. Note that the
-only particles that have hydrodynamical forces calculated between them are those
-in ``PartType0``.
+There are several groups that contain 'auxilliary' information, such as
+``Header``.  Particle data is placed in separate groups depending of the type of
+the particles. Some types are currently ignored by SWIFT but are kept in the
+file format for compatibility reasons.
+
++---------------------+------------------------+----------------------------+
+| HDF5 Group Name     | Physical Particle Type | In code ``enum part_type`` |
++=====================+========================+============================+
+| ``/PartType0/``     | Gas                    | ``swift_type_gas``         |
++---------------------+------------------------+----------------------------+
+| ``/PartType1/``     | Dark Matter            | ``swift_type_dark_matter`` |
++---------------------+------------------------+----------------------------+
+| ``/PartType2/``     | Ignored                |                            |
++---------------------+------------------------+----------------------------+
+| ``/PartType3/``     | Ignored                |                            |
++---------------------+------------------------+----------------------------+
+| ``/PartType4/``     | Stars                  | ``swift_type_star``        |
++---------------------+------------------------+----------------------------+
+| ``/PartType5/``     | Black Holes            | ``swift_type_black_hole``  |
++---------------------+------------------------+----------------------------+
+
+The last column in the table gives the ``enum`` value from ``part_type.h``
+corresponding to a given entry in the files.
+
+Note that the only particles that have hydrodynamical forces calculated between
+them are those in ``PartType0``.
 
 
 Necessary Components
@@ -55,27 +73,34 @@ script.
 Header
 ~~~~~~
 
-In ``Header``, the following attributes are required:
-
-+ ``BoxSize``, a floating point number or N-dimensional (usually 3) array
-  that describes the size of the box.
+In the ``/Header/`` group, the following attributes are required:
+
++ ``Dimension``, an integer indicating the dimensionality of the ICs (1,2 or 3).
+  Note that this parameter is an addition to the GADGET-2 format and will be
+  ignored by GADGET. SWIFT will use this value to verify that the dimensionality
+  of the code matches the ICs. If this parameter is not provided, it defaults
+  to 3.
++ ``BoxSize``, a floating point number or N-dimensional (usually 3) array that
+  describes the size of the box. If only one number is provided (as per the
+  GADGET-2 standard) then the box is assumed have the same size along all the
+  axis. In cosmological runs, this is the comoving box-size expressed in the
+  units specified in the ``/Units`` group (see below). Note that, unlike GADGET,
+  we express all quantities in "h-free" units. So that, for instance, we express
+  the box side-length in ``Mpc`` and not ``Mpc/h``. 
++ ``NumPart_Total``, a length 6 array of integers that tells the code how many
+  particles of each type are in the initial conditions file. Unlike traditional
+  GADGET-2 files, these can be >2^31.
++ ``NumPart_Total_HighWord``, a historical length-6 array that tells the code
+  the number of high word particles in the initial conditions there are. If you
+  are unsure, just set this to ``[0, 0, 0, 0, 0, 0]``. This does have to be
+  present but can be a set of 0s unless you have more than 2^31 particles and
+  want to be fully compliant with GADGET-2. Note that, as SWIFT supports
+  ``NumPart_Total`` to be >2^31, the use of ``NumPart_Total_HighWord`` is only
+  here for compatibility reasons.
 + ``Flag_Entropy_ICs``, a historical value that tells the code if you have
   included entropy or internal energy values in your intial conditions files.
-  Acceptable values are 0 or 1.
-+ ``NumPart_Total``, a length 6 array of integers that tells the code how many
-  particles are of each type are in the initial conditions file.
-+ ``NumPart_Total_HighWord``, a historical length-6 array that tells the code 
-  the number of high word particles in the initial conditions there are. If
-  you are unsure, just set this to ``[0, 0, 0, 0, 0, 0]``. This does have to be
-  present, but unlike GADGET-2, this can be a set of 0s unless you have more than
-  2^31 particles.
-+ ``NumFilesPerSnapshot``, again a historical integer value that tells the code
-  how many files there are per snapshot. You will probably want to set this to 1
-  and simply have a single HDF5 file for your initial conditions; SWIFT can
-  leverage parallel-HDF5 to read from this single file in parallel.
-+ ``NumPart_ThisFile``, a length 6 array of integers describing the number of
-  particles in this file. If you have followed the above advice, this will be
-  exactly the same as the ``NumPart_Total`` array.
+  Acceptable values are 0 or 1. We recommend using internal energies over
+  entropy in the ICs and hence have this flag set to 0.
 
 You may want to include the following for backwards-compatibility with many
 GADGET-2 based analysis programs:
@@ -83,63 +108,102 @@ GADGET-2 based analysis programs:
 + ``MassTable``, an array of length 6 which gives the masses of each particle
   type. SWIFT ignores this and uses the individual particle masses, but some
   programs will crash if it is not included.
-+ ``Time``, the internal code time of the start (set this to 0).
-
++ ``NumPart_ThisFile``, a length 6 array of integers describing the number of
+  particles in this file. If you have followed the above advice, this will be
+  exactly the same as the ``NumPart_Total`` array. As SWIFT only uses ICs
+  contained in a single file, this is not necessary for SWIFT-only ICs.
++ ``NumFilesPerSnapshot``, again a historical integer value that tells the code
+  how many files there are per snapshot. You will probably want to set this to 1.
++ ``Time``, time of the start of the simulation in internal units or expressed
+  as a scale-factor for cosmological runs. SWIFT ignores this and reads it from
+  the parameter file.
+  
 RuntimePars
 ~~~~~~~~~~~
 
-In ``RuntimePars``, the following attributes are required:
+In the ``/RuntimePars/``, the following attributes are required:
 
 + ``PeriodicBoundaryConditionsOn``, a flag to tell the code whether or not you
   have periodic boundaries switched on. Again, this is historical; it should be
   set to 1 (default) if you have the code running in periodic mode, or 0 otherwise.
 
 
-Units
-~~~~~
-
-In ``Units``, you will need to specify what units your initial conditions are
-in. If these are not present, the code assumes that you are using the same
-units for your initial conditions as are in your parameterfile, but it is best
-to include them to be on the safe side. You will need:
-
-+ ``Unit current in cgs (U_I)``
-+ ``Unit length in cgs (U_L)``
-+ ``Unit mass in cgs (U_M)``
-+ ``Unit temperature in cgs (U_T)``
-+ ``Unit time in cgs (U_t)``
-
-These are all floating point numbers.
-
-
 Particle Data
 ~~~~~~~~~~~~~
 
 Now for the interesting part! You can include particle data groups for each
-individual particle type (e.g. ``PartType0``) that have the following _datasets_:
+individual particle type (e.g. ``/PartType0/``) that have the following *datasets*:
 
 + ``Coordinates``, an array of shape (N, 3) where N is the number of particles
-  of that type, that are the cartesian co-ordinates of the particles. Co-ordinates
-  must be positive, but will be wrapped on reading to be within the periodic box.
-+ ``Velocities``, an array of shape (N, 3) that is the cartesian velocities 
-  of the particles.
+  of that type, that are the cartesian co-ordinates of the
+  particles. Co-ordinates must be within the box so, in the case of a cube
+  within [0, L)^3 where L is the side-length of the simulation volume. In the
+  case of cosmological simulations, these are the co-moving positions.
++ ``Velocities``, an array of shape (N, 3) that is the cartesian velocities of
+  the particles. When running cosmological simulations, these are the peculiar
+  velocities. Note that this is different from GADGET which uses peculiar
+  velocities divided by ``sqrt(a)`` (see below for a fix).
 + ``ParticleIDs``, an array of length N that are unique identifying numbers for
   each particle. Note that these have to be unique to a particle, and cannot be
-  the same even between particle types. Please ensure that your IDs are positive
-  integer numbers.
+  the same even between particle types. The **IDs must be >1**. 0 or negative
+  IDs will be rejected by the code.
 + ``Masses``, an array of length N that gives the masses of the particles.
 
-For ``PartType0`` (i.e. particles that interact through hydrodynamics), you will
+For ``PartType0`` (i.e. particles that interact through hydro-dynamics), you will
 need the following auxilliary items:
 
-+ ``InternalEnergy``, an array of length N that gives the internal energies of
-  the particles. For PressureEntropy, you can specify ``Entropy`` instead.
-+ ``SmoothingLength``, the smoothing lenghts of the particles. These will be
-  tidied up a bit, but it is best if you provide accurate numbers.
++ ``SmoothingLength``, the smoothing lengths of the particles. These will be
+  tidied up a bit, but it is best if you provide accurate numbers. In
+  cosmological runs, these are the co-moving smoothing lengths.
++ ``InternalEnergy``, an array of length N that gives the internal energies per
+  unit mass of the particles. If the hydro-scheme used in the code is based on
+  another thermodynamical quantity (entropy or total energy, etc.), the
+  conversion will happen inside the code. In cosmological runs, this is the
+  **physical** internal energy per unit mass. This has the dimension of velocity
+  squared.
+
+  
+Note that for cosmological runs, all quantities have to be expressed in "h-free"
+dimensions. This means ``Mpc`` and not ``Mpc/h`` for instance. If the ICs have
+been generated for GADGET (where h-full values are expected), the parameter
+``InitialConditions:cleanup_h_factors`` can be set to ``1`` in the
+:ref:`Parameter_File_label` to make SWIFT convert the quantities read in to
+h-free quantities. Switching this parameter on will also affect the box size
+read from the ``/Header/`` group (see above).
+
+Similarly, GADGET cosmological ICs have traditionally used velocities expressed
+as peculiar velocities divided by ``sqrt(a)``. This can be undone by swicthing
+on the parameter ``InitialConditions:cleanup_velocity_factors`` in the
+:ref:`Parameter_File_label`.
+
+     
+Optional Components
+-------------------
+
+In the ``/Units/`` HDF5 group, you cans specify what units your initial conditions are
+in. If this group is not present, the code assumes that you are using the same
+units for your initial conditions as in your :ref:`Parameter_File_label`
+(i.e. as the internal units system used by the code), but it is best to include
+them to be on the safe side. You will need:
+
++ ``Unit length in cgs (U_L)``
++ ``Unit mass in cgs (U_M)``
++ ``Unit time in cgs (U_t)``
++ ``Unit current in cgs (U_I)``
++ ``Unit temperature in cgs (U_T)``
+
+These are all floating point numbers. Note that we specify the time units and
+not the velocity units.
+
+If the units specified in the initial conditions are different from the internal
+units (specified in the parameter file), SWIFT will perform a conversion of all
+the quantities when reading in the ICs. This includes a conversion of the box
+size read from the ``/Header/`` group.
 
 
+     
 Summary
-~~~~~~~
+-------
 
 You should have an HDF5 file with the following structure:
 
@@ -147,11 +211,9 @@ You should have an HDF5 file with the following structure:
 
    Header/
      BoxSize=[x, y, z]
-     Flag_Entropy_ICs=1
-     NumPart_Total=[0, 1, 2, 3, 4, 5]
+     Flag_Entropy_ICs=0
+     NumPart_Total=[0, 1, 0, 0, 4, 5]
      NumPart_Total_HighWord=[0, 0, 0, 0, 0, 0]
-     NumFilesPerSnapshot=1
-     NumPart_ThisFile=[0, 1, 2, 3, 4, 5]
    RuntimePars/
      PeriodicBoundariesOn=1
    Units/

examples/AgoraDisk/agora_disk.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e21 # kpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
+
+Scheduler:
+  max_top_level_cells: 8
+  
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   0.5  # The end time of the simulation (in internal units).
+  dt_min:     1e-10 # The minimal 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
+Snapshots:
+  basename:            agora_disk # Common part of the name of output files
+  time_first:          0.    # Time of the first output (in internal units)
+  delta_time:          1e-2  # Time difference between consecutive outputs (in internal units)
+  compression:         4
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-3 # Time between statistics output
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                    0.025    # Constant dimensionless multiplier for time integration.
+  theta:                  0.7     # Opening angle (Multipole acceptance criterion)
+  comoving_softening:     0.08 # Comoving softening length (in internal units).
+  max_physical_softening: 0.08    # Physical softening length (in internal units).
+  mesh_side_length:       32        # Number of cells along each axis for the periodic gravity mesh.
+  
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10      # (internal units)
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./agora_disk.hdf5     # The file to read
+  cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
+  shift:    [674.1175, 674.1175, 674.1175]   # (Optional) A shift to apply to all particles read from the ICs (in internal units).
+
+# Dimensionless pre-factor for the time-step condition
+LambdaCooling:
+  lambda_cgs:                  1.0e-22    # Cooling rate (in cgs units)
+  minimum_temperature:         1.0e2      # Minimal temperature (Kelvin)
+  mean_molecular_weight:       0.59       # Mean molecular weight
+  hydrogen_mass_abundance:     0.75       # Hydrogen mass abundance (dimensionless)
+  cooling_tstep_mult:          1.0        # Dimensionless pre-factor for the time-step condition
+
+# Cooling with Grackle 2.0
+GrackleCooling:
+  CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+  WithUVbackground: 1 # Enable or not the UV background
+  Redshift: 0 # Redshift to use (-1 means time based redshift)
+  WithMetalCooling: 1 # Enable or not the metal cooling
+  ProvideVolumetricHeatingRates: 0 # User provide volumetric heating rates
+  ProvideSpecificHeatingRates: 0 # User provide specific heating rates
+  SelfShieldingMethod: 0 # Grackle (<= 3) or Gear self shielding method
+  OutputMode: 1 # Write in output corresponding primordial chemistry mode
+  MaxSteps: 1000
+  ConvergenceLimit: 1e-2

examples/AgoraDisk/changeType.py

+#!/usr/bin/env python3
+
+from h5py import File
+from sys import argv
+import numpy as np
+
+"""
+Change particle types in order to match the implemented types
+"""
+
+# number of particle type
+N_type = 6
+
+debug = 0
+
+
+def getOption():
+    if len(argv) != 2:
+        raise IOError("You need to provide a filename")
+
+    # get filename and read it
+    filename = argv[-1]
+
+    return filename
+
+
+def groupName(part_type):
+    return "PartType%i" % part_type
+
+
+def changeType(f, old, new):
+    # check if directory exists
+    old_group = groupName(old)
+    if old_group not in f:
+        raise IOError("Cannot find group '%s'" % old)
+    old = f[old_group]
+
+    new_group = groupName(new)
+    if new_group not in f:
+        f.create_group(new_group)
+    new = f[new_group]
+
+    for name in old:
+        if debug:
+            print("Moving '%s' from '%s' to '%s'"
+                  % (name, old_group, new_group))
+
+        tmp = old[name][:]
+        del old[name]
+        if name in new:
+            new_tmp = new[name][:]
+            if debug:
+                print("Found previous data:", tmp.shape, new_tmp.shape)
+            tmp = np.append(tmp, new_tmp, axis=0)
+            del new[name]
+
+        if debug:
+            print("With new shape:", tmp.shape)
+        new.create_dataset(name, tmp.shape)
+        new[name][:] = tmp
+
+    del f[old_group]
+
+
+def countPart(f):
+    npart = []
+    for i in range(N_type):
+        name = groupName(i)
+        if name in f:
+            grp = f[groupName(i)]
+            N = grp["Masses"].shape[0]
+        else:
+            N = 0
+        npart.append(N)
+
+    f["Header"].attrs["NumPart_ThisFile"] = npart
+    f["Header"].attrs["NumPart_Total"] = npart
+    f["Header"].attrs["NumPart_Total_HighWord"] = [0]*N_type
+
+
+if __name__ == "__main__":
+    filename = getOption()
+
+    f = File(filename)
+
+    changeType(f, 2, 1)
+    changeType(f, 3, 1)
+    changeType(f, 4, 1)
+
+    countPart(f)
+
+    f.close()

examples/AgoraDisk/cleanupSwift.py

+#!/usr/bin/env python3
+
+# ./translate_particles.py filename output_name
+from h5py import File
+import sys
+from shutil import copyfile
+
+NPartType = 1
+filename = sys.argv[-2]
+out = sys.argv[-1]
+
+copyfile(filename, out)
+
+f = File(out)
+
+name = "PartType0/ElementAbundance"
+if name in f:
+    del f[name]
+
+for i in range(NPartType):
+    name = "PartType%i" % i
+    if name not in f:
+        continue
+
+    grp = f[name + "/SmoothingLength"]
+    grp[:] *= 1.823
+
+f.close()

examples/AgoraDisk/getIC.sh

+#!/bin/bash
+
+if [ "$#" -ne 1 ]; then
+    echo "You need to provide the resolution (e.g. ./getIC.sh low)."
+    echo "The possible options are low, med and high."
+    exit
+fi
+
+wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/$1

examples/AgoraDisk/getSolution.sh

+#!/bin/bash
+
+OPTIND=1
+
+with_cooling=0
+
+function show_help {
+    echo "Valid options are:"
+    echo "\t -h \t Show this help"
+    echo "\t -C \t Download solution with cooling"
+}
+
+while getopts "h?C" opt; do
+    case "$opt" in
+	h|\?)
+	    show_help
+	    exit
+	    ;;
+	C)
+	    with_cooling=1
+	    ;;
+    esac
+done
+
+# cleanup work space
+rm snapshot_0000.hdf5 snapshot_0500.hdf5
+
+if [ $with_cooling -eq 1 ]; then
+    wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/with_cooling/snapshot_0000.hdf5
+    wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/with_cooling/snapshot_0500.hdf5
+else
+    wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/without_cooling/snapshot_0000.hdf5
+    wget https://obswww.unige.ch/~lhausamm/swift/IC/AgoraDisk/Gear/without_cooling/snapshot_0500.hdf5
+fi
+
+

examples/AgoraDisk/run.sh

+#!/bin/bash
+
+# This example is based on the AGORA disk article (DOI: 10.3847/1538-4357/833/2/202)
+
+# currently only the low resolution is available
+sim=low
+
+# enable cooling or not
+cooling=0
+
+# make run.sh fail if a subcommand fails
+set -e
+
+# define flags
+flag=
+if [ $cooling -eq 1 ]
+then
+    flag=-C
+fi
+
+# Generate the initial conditions if they are not present.
+if [ ! -e $sim.hdf5 ]
+then
+    echo "Fetching initial glass file for the Sedov blast example..."
+    ./getIC.sh $sim.hdf5
+fi
+
+# Get the Grackle cooling table
+if [ ! -e CloudyData_UVB=HM2012.h5 ]
+then
+    echo "Fetching the Cloudy tables required by Grackle..."
+    ../getCoolingTable.sh
+fi
+
+# copy the initial conditions
+cp $sim.hdf5 agora_disk.hdf5
+# Update the particle types
+python3 changeType.py agora_disk.hdf5
+
+# Run SWIFT
+#../swift $flag -s -G -t 4 agora_disk.yml 2>&1 | tee output.log
+
+
+echo "Changing smoothing length to be Gadget compatible"
+python3 cleanupSwift.py agora_disk_0000.hdf5 agora_disk_IC.hdf5
+python3 cleanupSwift.py agora_disk_0050.hdf5 agora_disk_500Myr.hdf5
+
+echo "Fetching GEAR solution..."
+./getSolution.sh $flag
+
+echo "Plotting..."
+python3 plotSolution.py $flag

examples/BigCosmoVolume/makeIC.py

-###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # 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 urllib
-import os.path
-import numpy as np
-import sys
-
-outputName = "bigCosmoVolume.hdf5" 
-
-
-#--------------------------------------------------
-if len(sys.argv) != 3:
-    print "Invalid number of arguments. Need to provide down-sampling factor [0.,1.] and number of copies (integer)"
-    print "Example: python makeIC.py 0.8 4"
-    exit()
-
-downsample = float(sys.argv[1])
-n_copy = int(sys.argv[2])
-
-if n_copy < 1:
-    print "Number of copy must be >1"
-    exit()
-
-if downsample > 1. or downsample <= 0.:
-    print "Down-sampling factor must be in [0,1]"
-    exit()
-
-#--------------------------------------------------
-
-# Download the tile
-if (not os.path.isfile("tile.hdf5")):
-    print "Downloading initial tile..."
-    urllib.urlretrieve ("http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/tile.hdf5", "tile.hdf5")
-    print "Done."
-else:
-    print "Tile already exists. No need to download..."
-
-# Read in the tile
-inputFile = h5py.File("tile.hdf5", 'r+')
-
-grp = inputFile["/Header"]
-boxSize = grp.attrs["BoxSize"]
-numPart = grp.attrs["NumPart_Total"][0]
-
-coords = inputFile["/PartType0/Coordinates"][:,:]
-v = inputFile["/PartType0/Velocities"][:,:]
-m = inputFile["/PartType0/Masses"][:]
-h = inputFile["/PartType0/SmoothingLength"][:]
-u = inputFile["/PartType0/InternalEnergy"][:]
-ids = np.array(range(np.size(u)), dtype='L') + 1
-
-# Downsample
-print "Downsampling..."
-indices = np.array(range(np.size(ids)))
-np.random.shuffle(indices)
-
-numPart *= downsample
-indices = indices < numPart
-
-coords = coords[indices,:]
-v = v[indices,:]
-m = m[indices]
-h = h[indices] / 1.825742 # Correct from Gadget defintion of h to physical definition
-u = u[indices]
-ids = ids[indices]
-
-numPart = np.size(ids)
-
-# Now replicate the tile
-if n_copy > 1:
-
-    print "Tiling..."
-
-    coords_tile = np.copy(coords)
-    v_tile = np.copy(v)
-    m_tile = np.copy(m)
-    h_tile = np.copy(h)
-    u_tile = np.copy(u)
-    ids_tile = np.copy(ids)
-
-    coords = np.zeros((0,3))
-    v = np.zeros((0,3))
-    m = np.zeros(0)
-    h = np.zeros(0)
-    u = np.zeros(0)
-    ids = np.zeros(0, dtype='L')
-
-    count = 0
-
-    for i in range(n_copy):
-        for j in range(n_copy):
-            for k in range(n_copy):
-                coords = np.append(coords, coords_tile + np.array([ i * boxSize[0], j * boxSize[1], k * boxSize[2] ]), axis=0)
-                v = np.append(v, v_tile, axis=0)
-                m = np.append(m, m_tile)
-                h = np.append(h, h_tile)
-                u = np.append(u, u_tile)
-                
-                ids = np.append(ids, ids_tile + count*numPart)
-                count+=1
-
-
-    numPart *= n_copy**3
-    boxSize *= n_copy
-
-# Copy the tile out
-file = h5py.File(outputName, 'w')
-
-# Header
-grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 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
-grp.attrs["Dimension"] = 3
-
-#Runtime parameters
-grp = file.create_group("/RuntimePars")
-grp.attrs["PeriodicBoundariesOn"] = 1
-
-#Particle group
-grp = file.create_group("/PartType0")
-grp.create_dataset('Coordinates', data=coords, dtype='d')
-grp.create_dataset('Velocities', data=v, dtype='f')
-grp.create_dataset('Masses', data=m, dtype='f')
-grp.create_dataset('SmoothingLength', data=h, dtype='f')
-grp.create_dataset('InternalEnergy', data=u, dtype='f')
-grp.create_dataset('ParticleIDs', data=ids, dtype='L')
-
-
-

examples/BigPerturbedBox/makeIC_fcc.py

-###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # 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 random
-from numpy import *
-
-# Generates a swift IC file for the Sedov blast test in a periodic cubic box
-
-# Parameters
-periodic= 1      # 1 For periodic box
-boxSize = 10.
-L = 128           # Number of particles boxes along one axis
-rho = 1.          # Density
-P = 1.e-5         # Pressure
-E0= 1.e2          # Energy of the explosion
-pert = 0.025
-gamma = 5./3.     # Gas adiabatic index
-fileName = "perturbedBox.hdf5" 
-
-
-#---------------------------------------------------
-numPart = 4*(L**3) 
-mass = boxSize**3 * rho / numPart
-internalEnergy = P / ((gamma - 1.)*rho)
-off = array( [ [ 0.0 , 0.0 , 0.0 ] , [ 0.0 , 0.5 , 0.5 ] , [ 0.5 , 0.0 , 0.5 ] , [ 0.5 , 0.5 , 0.0 ] ] );
-hbox = boxSize / L
-
-# if L%2 == 0:
-#     print "Number of particles along each dimension must be odd."
-#     exit()
-
-#Generate particles
-coords = zeros((numPart, 3))
-v      = zeros((numPart, 3))
-m      = zeros((numPart, 1))
-h      = zeros((numPart, 1))
-u      = zeros((numPart, 1))
-ids    = zeros((numPart, 1), dtype='L')
-
-for i in range(L):
-    for j in range(L):
-        for k in range(L):
-            x = (i + 0.25) * hbox
-            y = (j + 0.25) * hbox
-            z = (k + 0.25) * hbox
-            for ell in range(4): 
-                index = 4*(i*L*L + j*L + k) + ell
-                coords[index,0] = x + off[ell,0] * hbox
-                coords[index,1] = y + off[ell,1] * hbox
-                coords[index,2] = z + off[ell,2] * hbox
-                v[index,0] = 0.
-                v[index,1] = 0.
-                v[index,2] = 0.
-                m[index] = mass
-                h[index] = 2.251 / 4 * hbox
-                u[index] = internalEnergy
-                ids[index] = index
-                coords[index,0] += random.random() * pert * hbox
-                coords[index,1] += random.random() * pert * hbox
-                coords[index,2] += random.random() * pert * hbox
-
-
-#--------------------------------------------------
-
-#File
-file = h5py.File(fileName, 'w')
-
-# Header
-grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] =  [numPart, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 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
-
-#Runtime parameters
-grp = file.create_group("/RuntimePars")
-grp.attrs["PeriodicBoundariesOn"] = periodic
-
-#Particle group
-grp = file.create_group("/PartType0")
-ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
-ds[()] = coords
-ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
-ds[()] = v
-ds = grp.create_dataset('Masses', (numPart,1), 'f')
-ds[()] = m
-ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
-ds[()] = h
-ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
-ds[()] = u
-ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
-ds[()] = ids + 1
-
-file.close()