Merge branch 'master' into root-level-parts

 Conflicts:
	src/partition.h

Parameter handling imported from old main.c.

.gitignore

@@ -25,17 +25,30 @@ examples/swift_mindt
 examples/swift_mindt_mpi
 examples/swift_mpi
 
-tests/testVectorize
-tests/brute_force.dat
-tests/swift_dopair.dat
+tests/testPair
+tests/brute_force_standard.dat
+tests/swift_dopair_standard.dat
+tests/brute_force_perturbed.dat
+tests/swift_dopair_perturbed.dat
+tests/test27cells
+tests/brute_force_27_standard.dat
+tests/swift_dopair_27_standard.dat
+tests/brute_force_27_perturbed.dat
+tests/swift_dopair_27_perturbed.dat
 tests/testGreetings
 tests/testReading
 tests/input.hdf5
 tests/testSingle
 tests/testTimeIntegration
 tests/testSPHStep
+tests/testKernel
+tests/testParser
+tests/parser_output.yml
 
 theory/latex/swift.pdf
+theory/kernel/kernels.pdf
+theory/kernel/kernel_derivatives.pdf
+theory/kernel/kernel_definitions.pdf
 
 m4/libtool.m4
 m4/ltoptions.m4

README

- Welcome to the cosmological code
+ Welcome to the cosmological hydrodynamical code
     ______       _________________
    / ___/ |     / /  _/ ___/_  __/
    \__ \| | /| / // // /_   / /   
@@ -6,8 +6,26 @@
  /____/ |__/|__/___/_/    /_/     
  SPH With Inter-dependent Fine-grained Tasking
 
-Website: www.swiftsim.com
-Twitter: @SwiftSimulation
+ Website: www.swiftsim.com
+ Twitter: @SwiftSimulation
 
-See INSTALL.swift for instructions.
+See INSTALL.swift for install instructions.
 
+Usage: swift [OPTION] PARAMFILE
+
+Valid options are:
+  -c          Run with cosmological time integration
+  -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.
+  -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
+  -s          Run with SPH
+  -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
+
+See the file examples/parameter_example.yml for an example of parameter file.

configure.ac

@@ -287,11 +287,11 @@ AC_SUBST([METIS_LIBS])
 AC_SUBST([METIS_INCS])
 AM_CONDITIONAL([HAVEMETIS],[test -n "$METIS_LIBS"])
 
-# Check for zlib.
-AC_CHECK_LIB([z],[gzopen],[
-    AC_DEFINE([HAVE_LIBZ],[1],[Set to 1 if zlib is installed.])
-    LDFLAGS="$LDFLAGS -lz"
-    ],[])
+# # Check for zlib.
+# AC_CHECK_LIB([z],[gzopen],[
+#     AC_DEFINE([HAVE_LIBZ],[1],[Set to 1 if zlib is installed.])
+#     LDFLAGS="$LDFLAGS -lz"
+#     ],[])
 
 
 # Check for HDF5. This is required.

doc/Doxyfile.in

@@ -759,7 +759,9 @@ WARN_LOGFILE           =
 # spaces.
 # Note: If this tag is empty the current directory is searched.
 
-INPUT                  = @top_srcdir@ @top_srcdir@/src @top_srcdir@/src/hydro/Minimal @top_srcdir@/src/gravity/Default
+INPUT                  =  @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
+INPUT		       += @top_srcdir@/src/hydro/Minimal @top_srcdir@/src/gravity/Default
+INPUT		       += @top_srcdir@/src/riemann 
 
 # 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

examples/BigCosmoVolume/makeIC.py

@@ -77,7 +77,7 @@ indices = indices < numPart
 coords = coords[indices,:]
 v = v[indices,:]
 m = m[indices]
-h = h[indices]
+h = h[indices] / 1.825742 # Correct from Gadget defintion of h to physical definition
 u = u[indices]
 ids = ids[indices]
 

examples/CosmoVolume/cosmoVolume.yml

+
+# Define the system of units to use internally. 
+UnitSystem:
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters for the task scheduling
+Scheduler:
+  nr_threads:       16        # The number of threads per MPI rank to use.
+  nr_queues:        0        # The number of task queues to use. Use 0  to let the system decide.
+  cell_max_size:    8000000  # Maximal number of interactions per task (this is the default value).
+  cell_sub_size:    5000     # Maximal number of interactions per sub-task  (this is the default value).
+  cell_split_size:  400      # Maximal number of particles per cell (this is the default value).
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.    # 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_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      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:  0.6      # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./cosmoVolume.hdf5     # The file to read
+  h_scaling:  1.                    # A scaling factor to apply to all smoothing lengths in the ICs.
+  shift_x:    0.                    # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    0.
+  shift_z:    0.
+
+# Parameters govering domain decomposition
+DomainDecomposition:
+  initial_type:       m     # The initial strategy ("g", "m", "w", or "v"). See documentation for details.
+  initial_grid_x:    10     # Grid size if the 'g' strategy is chosen.
+  initial_grid_y:    10
+  initial_grid_z:    10
+  repartition_type:   b     # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
+ 

examples/CosmoVolume/run.sh

+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e cosmoVolume.hdf5 ]
+then
+    echo "Fetching initial conditions for the cosmo volume example..."
+    ./getIC.sh
+fi
+
+../swift -s cosmoVolume.yml

examples/GreshoVortex/makeIC.py

@@ -30,6 +30,7 @@ factor = 3
 boxSize = [ 1.0 , 1.0, 1.0/factor ]
 L = 120           # Number of particles along one axis
 gamma = 5./3.     # Gas adiabatic index
+eta = 1.2349      # 48 ngbs with cubic spline kernel
 rho = 1           # Gas density
 P0 = 0.           # Constant additional pressure (should have no impact on the dynamics)
 fileName = "greshoVortex.hdf5" 
@@ -73,7 +74,7 @@ for i in range(L):
             v[index,1] =  v_phi * (x - boxSize[0] / 2) / r
             v[index,2] = 0.
             m[index] = mass
-            h[index] = 2.251 * boxSize[0] / L
+            h[index] = eta * boxSize[0] / L
             P = P0
             if r < 0.2:
                 P = P + 5. + 12.5*r2
@@ -105,6 +106,14 @@ grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 #Particle group
 grp = file.create_group("/PartType0")
 ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')

examples/MultiTypes/makeIC.py

@@ -32,6 +32,7 @@ Lgas = int(sys.argv[1])  # Number of particles along one axis
 rhoGas = 2.              # Density
 P = 1.                   # Pressure
 gamma = 5./3.            # Gas adiabatic index
+eta = 1.2349             # 48 ngbs with cubic spline kernel
 rhoDM = 1.
 Ldm = int(sys.argv[2])  # Number of particles along one axis
 
@@ -61,11 +62,18 @@ grp.attrs["NumFilesPerSnapshot"] = 1
 grp.attrs["MassTable"] = [0.0, massDM, 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
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 
 # Gas Particle group
 grp = file.create_group("/PartType0")
@@ -80,7 +88,7 @@ ds = grp.create_dataset('Masses', (numGas,1), 'f')
 ds[()] = m
 m = zeros(1)
 
-h = full((numGas, 1), 1.1255 * boxSize / Lgas)
+h = full((numGas, 1), eta * boxSize / Lgas)
 ds = grp.create_dataset('SmoothingLength', (numGas,1), 'f')
 ds[()] = h
 h = zeros(1)

examples/PerturbedBox/makeIC.py

@@ -90,6 +90,14 @@ grp.attrs["NumPart_Total"] = numPart
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 #Particle group
 grp = file.create_group("/PartType0")
 ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')

examples/SedovBlast/makeIC.py

@@ -33,6 +33,7 @@ P = 1.e-5         # Pressure
 E0= 1.e2          # Energy of the explosion
 pert = 0.1
 gamma = 5./3.     # Gas adiabatic index
+eta = 1.2349      # 48 ngbs with cubic spline kernel
 fileName = "sedov.hdf5" 
 
 
@@ -67,7 +68,7 @@ for i in range(L):
             v[index,1] = 0.
             v[index,2] = 0.
             m[index] = mass
-            h[index] = 1.1255 * boxSize / L
+            h[index] = eta * boxSize / L
             u[index] = internalEnergy
             ids[index] = index + 1
             if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 2.01 * boxSize/L:
@@ -98,6 +99,14 @@ grp.attrs["Flag_Entropy_ICs"] = 0
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 #Particle group
 grp = file.create_group("/PartType0")
 grp.create_dataset('Coordinates', data=coords, dtype='d')

examples/SedovBlast/makeIC_fcc.py

@@ -33,6 +33,7 @@ P = 1.e-5         # Pressure
 E0= 1.e2          # Energy of the explosion
 pert = 0.025
 gamma = 5./3.     # Gas adiabatic index
+eta = 1.2349          # 48 ngbs with cubic spline kernel
 fileName = "sedov.hdf5" 
 
 
@@ -70,7 +71,7 @@ for i in range(L):
                 v[index,1] = 0.
                 v[index,2] = 0.
                 m[index] = mass
-                h[index] = 1.1255 * hbox
+                h[index] = eta * hbox
                 u[index] = internalEnergy
                 ids[index] = index + 1
                 if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 1.2 * hbox:
@@ -101,6 +102,14 @@ grp.attrs["Flag_Entropy_ICs"] = 0
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 #Particle group
 grp = file.create_group("/PartType0")
 grp.create_dataset('Coordinates', data=coords, dtype='d')

examples/SedovBlast/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e sedov.hdf5 ]
+then
+    echo "Generating initial conditions for the SedovBlast example..."
+    python makeIC_fcc.py
+fi
+
+../swift -s sedov.yml

examples/SedovBlast/sedov.yml

+
+# Define the system of units to use internally. 
+UnitSystem:
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters for the task scheduling
+Scheduler:
+  nr_threads:       16       # The number of threads per MPI rank to use.
+  nr_queues:        0        # The number of task queues to use. Use 0  to let the system decide.
+  cell_max_size:    8000000  # Maximal number of interactions per task (this is the default value).
+  cell_sub_size:    5000     # Maximal number of interactions per sub-task  (this is the default value).
+  cell_split_size:  400      # Maximal number of particles per cell (this is the default value).
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.    # The end time of the simulation (in internal units).
+  dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      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:  1.       # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./sedov.hdf5          # The file to read
+  h_scaling:  1.                    # A scaling factor to apply to all smoothing lengths in the ICs.
+  shift_x:    0.                    # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    0.
+  shift_z:    0.
+
+# Parameters govering domain decomposition
+DomainDecomposition:
+  initial_type:       m     # The initial strategy ("g", "m", "w", or "v"). See documentation for details.
+  initial_grid_x:    10     # Grid size if the 'g' strategy is chosen.
+  initial_grid_y:    10
+  initial_grid_z:    10
+  repartition_type:   b     # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
+ 

examples/SodShock/makeIC.py

@@ -43,14 +43,14 @@ vol = boxSize[0] * boxSize[1] * boxSize[2]
 glass1 = h5py.File("glass_001.hdf5")
 pos1 = glass1["/PartType0/Coordinates"][:,:]
 pos1 = pos1 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
-
+glass_h1 = glass1["/PartType0/SmoothingLength"][:] / factor
 
 #Read in high density glass
 # glass2 = h5py.File("../Glass/glass_50000.hdf5")
 glass2 = h5py.File("glass_002.hdf5")
 pos2 = glass2["/PartType0/Coordinates"][:,:]
 pos2 = pos2 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
-
+glass_h2 = glass2["/PartType0/SmoothingLength"][:] / factor
 
 #Generate high density region
 rho1 = 1.
@@ -61,9 +61,10 @@ coord1 = append(coord1, coord1 + [0.25, 0, 0], 0)
 # coord1 = append(coord1, pos1 + [0, 0.5, 0.5], 0)
 N1 = size(coord1)/3
 v1 = zeros((N1, 3))
-h1 = ones(N1) * 2.251 * 0.5 * vol / (size(pos1)/3)**(1./3.)
 u1 = ones(N1) * P1 / ((gamma - 1.) * rho1)
 m1 = ones(N1) * vol * 0.5 * rho1 / N1
+h1 = append(glass_h1, glass_h1, 0)
+h1 = append(h1, h1, 0)
 
 #Generate low density region
 rho2 = 0.25
@@ -74,9 +75,10 @@ coord2 = append(coord2, coord2 + [0.25, 0, 0], 0)
 # coord2 = append(coord2, pos2 + [0, 0.5, 0.5], 0)
 N2 = size(coord2)/3
 v2 = zeros((N2, 3))
-h2 = ones(N2) * 2.251 * 0.5 * vol / (size(pos2)/3)**(1./3.)
 u2 = ones(N2) * P2 / ((gamma - 1.) * rho2)
 m2 = ones(N2) * vol * 0.5 * rho2 / N2
+h2 = append(glass_h2, glass_h2, 0)
+h2 = append(h2, h2, 0)
 
 #Merge arrays
 numPart = N1 + N2
@@ -89,8 +91,8 @@ ids = zeros(numPart, dtype='L')
 for i in range(1, numPart+1):
     ids[i-1] = i
 
-#Final operations
-h /= 2
+#Final operation since we come from Gadget-2 cubic spline ICs
+h /= 1.825752
 
 #File
 file = h5py.File(fileName, 'w')
@@ -110,6 +112,14 @@ grp.attrs["Flag_Entropy_ICs"] = 0
 grp = file.create_group("/RuntimePars")
 grp.attrs["PeriodicBoundariesOn"] = periodic
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 #Particle group
 grp = file.create_group("/PartType0")
 grp.create_dataset('Coordinates', data=coords, dtype='d')

examples/SodShock/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e sodShock.hdf5 ]
+then
+    echo "Generating initial conditions for the SodShock example..."
+    python makeIC.py
+fi
+
+../swift -s sodShock.yml

examples/SodShock/sodShock.yml

+
+# Define the system of units to use internally. 
+UnitSystem:
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters for the task scheduling
+Scheduler:
+  nr_threads:       16        # The number of threads per MPI rank to use.
+  nr_queues:        0        # The number of task queues to use. Use 0  to let the system decide.
+  cell_max_size:    8000000  # Maximal number of interactions per task (this is the default value).
+  cell_sub_size:    5000     # Maximal number of interactions per sub-task.
+  cell_split_size:  400      # Maximal number of particles per cell (this is the default value).
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.    # The end time of the simulation (in internal units).
+  dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      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:  0.01     # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./sodShock.hdf5       # The file to read
+  h_scaling:  1.                    # A scaling factor to apply to all smoothing lengths in the ICs.
+  shift_x:    0.                    # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    0.
+  shift_z:    0.
+
+# Parameters govering domain decomposition
+DomainDecomposition:
+  initial_type:       m     # The initial strategy ("g", "m", "w", or "v"). See documentation for details.
+  initial_grid_x:    10     # Grid size if the 'g' strategy is chosen.
+  initial_grid_y:    10
+  initial_grid_z:    10
+  repartition_type:   b     # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
+ 

examples/UniformBox/makeIC.py

@@ -32,6 +32,7 @@ L = int(sys.argv[1])  # Number of particles along one axis
 rho = 2.              # Density
 P = 1.                # Pressure
 gamma = 5./3.         # Gas adiabatic index
+eta = 1.2349          # 48 ngbs with cubic spline kernel
 fileName = "uniformBox.hdf5" 
 
 #---------------------------------------------------
@@ -55,11 +56,18 @@ 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
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
 #Particle group
 grp = file.create_group("/PartType0")
 
@@ -73,7 +81,7 @@ ds = grp.create_dataset('Masses', (numPart,1), 'f')
 ds[()] = m
 m = zeros(1)
 
-h = full((numPart, 1), 1.1255 * boxSize / L)
+h = full((numPart, 1), eta * boxSize / L)
 ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
 ds[()] = h
 h = zeros(1)

examples/UniformBox/makeICbig.py

@@ -32,6 +32,7 @@ N = int(sys.argv[2])  # Write N particles at a time to avoid requiring a lot of
 rho = 2.              # Density
 P = 1.                # Pressure
 gamma = 5./3.         # Gas adiabatic index
+eta = 1.2349      # 48 ngbs with cubic spline kernel
 fileName = "uniformBox_%d.hdf5"%L
 
 #---------------------------------------------------
@@ -62,11 +63,19 @@ 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
 
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+
 #Particle group
 grp = file.create_group("/PartType0")
 
@@ -89,7 +98,7 @@ for n in range(n_iterations):
     ds_m[offset:offset+N] = m
     m = zeros(1)
 
-    h = full((N, 1), 1.1255 * boxSize / L)
+    h = full((N, 1), eta * boxSize / L)
     ds_h[offset:offset+N] = h
     h = zeros(1)
 
@@ -122,7 +131,7 @@ m = full((remainder, 1), mass)
 ds_m[offset:offset+remainder] = m
 m = zeros(1)
 
-h = full((remainder, 1), 1.1255 * boxSize / L)
+h = full((remainder, 1), eta * boxSize / L)
 ds_h[offset:offset+remainder] = h
 h = zeros(1)
 
@@ -139,7 +148,7 @@ coords = zeros((remainder, 3))
 coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
 coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
 coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
-ds_x[offset:offset+remainder,:] = coords
+ods_x[offset:offset+remainder,:] = coords
 
 print "Done", offset+remainder,"/", numPart
 

examples/UniformBox/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e uniformBox.hdf5 ]
+then
+    echo "Generating initial conditions for the uniform box example..."
+    python makeIC.py 100
+fi
+
+../swift -s uniformBox.yml