Merge branch 'master' into mpi_and_ht_affinity

Conflicts:
        src/engine.c

.gitignore

@@ -11,7 +11,7 @@ config.sub
 ltmain.sh
 libtool
 
-src/version.h
+src/version_string.h
 swift*.tar.gz
 doc/doxyfile.stamp
 doc/html/
@@ -49,6 +49,7 @@ theory/latex/swift.pdf
 theory/kernel/kernels.pdf
 theory/kernel/kernel_derivatives.pdf
 theory/kernel/kernel_definitions.pdf
+theory/paper_pasc/pasc_paper.pdf
 
 m4/libtool.m4
 m4/ltoptions.m4

README

@@ -23,6 +23,7 @@ Valid options are:
   -g          Run with an external gravitational potential
   -G          Run with self-gravity
   -s          Run with SPH
+  -t    {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
   -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

configure.ac

@@ -16,7 +16,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 # Init the project.
-AC_INIT([SWIFT],[0.2.0])
+AC_INIT([SWIFT],[0.3.0])
 AC_CONFIG_SRCDIR([src/space.c])
 AC_CONFIG_AUX_DIR([.])
 AM_INIT_AUTOMAKE

examples/CosmoVolume/cosmoVolume.yml

 
 # Define the system of units to use internally. 
-UnitSystem:
+InternalUnitSystem:
   UnitMass_in_cgs:     1   # Grams
   UnitLength_in_cgs:   1   # Centimeters
   UnitVelocity_in_cgs: 1   # Centimeters per second
@@ -9,7 +9,6 @@ UnitSystem:
 
 # 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).
@@ -22,13 +21,25 @@ TimeIntegration:
   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 governing the snapshots
+Snapshots:
+  basename:            cosmo # Common part of the name of output files
+  time_first:          0.    # Time of the first output (in internal units)
+  delta_time:          0.05  # Time difference between consecutive outputs (in internal units)
+  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 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.
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
   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).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  max_volume_change:     2.       # Maximal allowed change of volume over one time-step
 
 # Parameters related to the initial conditions
 InitialConditions:

examples/CosmoVolume/run.sh

@@ -7,4 +7,4 @@ then
     ./getIC.sh
 fi
 
-../swift -s cosmoVolume.yml
+../swift -s -t 16 cosmoVolume.yml

examples/ExternalPointMass/externalPointMass.yml

+
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.9885e33     # Grams
+  UnitLength_in_cgs:   3.0856776e21  # Centimeters
+  UnitVelocity_in_cgs: 1e5           # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters for the task scheduling
+Scheduler:
+  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:    8000000  # 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-3  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            pointMass # Common part of the name of output files
+  time_first:          0.        # Time of the first output (in internal units)
+  delta_time:          0.02      # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1.9885e33     # Grams
+  UnitLength_in_cgs:   3.0856776e21  # Centimeters
+  UnitVelocity_in_cgs: 1e5           # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# 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).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
+  max_ghost_iterations:  30       # Maximal number of iterations allowed to converge towards the smoothing length.
+  max_smoothing_length:  10.      # Maximal smoothing length allowed (in internal units).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  max_volume_change:     2.       # Maximal allowed change of volume over one time-step
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  Sphere.hdf5           # The file to read
+  h_scaling:  1.                    # A scaling factor to apply to all smoothing lengths in the ICs.
+  shift_x:    50.                   # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    50.
+  shift_z:    50.
+
+# 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.
+ 
+
+# External potential parameters
+PointMass:
+  position_x:      50.     # location of external point mass in internal units
+  position_y:      50.
+  position_z:      50.	
+  mass:            1e10     # mass of external point mass in internal units
+

examples/ExternalPointMass/makeIC.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2016 John A. Regan (john.a.regan@durham.ac.uk)
+ #                    Tom Theuns (tom.theuns@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 sys
+import numpy
+import math
+import random
+
+# Generates a random distriution of particles, for motion in an external potnetial centred at (0,0,0)
+
+# physical constants in cgs
+NEWTON_GRAVITY_CGS  = 6.672e-8
+SOLAR_MASS_IN_CGS   = 1.9885e33
+PARSEC_IN_CGS       = 3.0856776e18
+
+# choice of units
+const_unit_length_in_cgs   =   (1000*PARSEC_IN_CGS)
+const_unit_mass_in_cgs     =   (SOLAR_MASS_IN_CGS)
+const_unit_velocity_in_cgs =   (1e5)
+
+print "UnitMass_in_cgs:     ", const_unit_mass_in_cgs 
+print "UnitLength_in_cgs:   ", const_unit_length_in_cgs
+print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
+
+# derived units
+const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
+const_G                = ((NEWTON_GRAVITY_CGS*const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/(const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs)))
+print 'G=', const_G
+
+
+# Parameters
+periodic= 1            # 1 For periodic box
+boxSize = 100.         # 
+Radius  = boxSize / 4. # maximum radius of particles
+G       = const_G 
+Mass    = 1e10         
+
+N       = int(sys.argv[1])  # Number of particles
+L       = N**(1./3.)
+
+# these are not used but necessary for I/O
+rho = 2.              # Density
+P = 1.                # Pressure
+gamma = 5./3.         # Gas adiabatic index
+fileName = "Sphere.hdf5" 
+
+
+#---------------------------------------------------
+numPart        = N
+mass           = 1
+internalEnergy = P / ((gamma - 1.)*rho)
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] =  [0, numPart, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [0, numPart, 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, 0, 0, 0, 0, 0]
+
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+#Particle group
+#grp0 = file.create_group("/PartType0")
+grp1 = file.create_group("/PartType1")
+#generate particle positions
+radius = Radius * (numpy.random.rand(N))**(1./3.) 
+ctheta = -1. + 2 * numpy.random.rand(N)
+stheta = numpy.sqrt(1.-ctheta**2)
+phi    =  2 * math.pi * numpy.random.rand(N)
+r      = numpy.zeros((numPart, 3))
+# r[:,0] = radius * stheta * numpy.cos(phi)
+# r[:,1] = radius * stheta * numpy.sin(phi)
+# r[:,2] = radius * ctheta
+r[:,0] = radius
+#
+speed  = numpy.sqrt(G * Mass / radius)
+v      = numpy.zeros((numPart, 3))
+omega  = speed / radius
+period = 2.*math.pi/omega
+print 'period = minimum = ',min(period), ' maximum = ',max(period)
+
+v[:,0] = -omega * r[:,1]
+v[:,1] =  omega * r[:,0]
+
+ds = grp1.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+v = numpy.zeros(1)
+
+m = numpy.full((numPart, ), mass)
+ds = grp1.create_dataset('Masses', (numPart,), 'f')
+ds[()] = m
+m = numpy.zeros(1)
+
+h = numpy.full((numPart, ), 1.1255 * boxSize / L)
+ds = grp1.create_dataset('SmoothingLength', (numPart,), 'f')
+ds[()] = h
+h = numpy.zeros(1)
+
+u = numpy.full((numPart, ), internalEnergy)
+ds = grp1.create_dataset('InternalEnergy', (numPart,), 'f')
+ds[()] = u
+u = numpy.zeros(1)
+
+
+ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False)
+ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
+ds[()] = ids
+
+ds = grp1.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = r
+
+file.close()

examples/ExternalPointMass/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e Sphere.hdf5 ]
+then
+    echo "Generating initial conditions for the point mass potential box example..."
+    python makeIC.py 10000
+fi
+
+../swift -g -t 2 externalPointMass.yml

examples/ExternalPointMass/test.pro

+;
+;  test energy / angular momentum conservation of test problem
+;
+@physunits
+
+indir    = '/gpfs/data/tt/Codes/Swift-git/swiftsim/examples/'
+basefile = 'output_'
+nfiles   = 657
+nfollow  = 100 ; number of particles to follow
+eout     = fltarr(nfollow, nfiles)
+ekin     = fltarr(nfollow, nfiles)
+epot     = fltarr(nfollow, nfiles)
+tout     = fltarr(nfiles)
+; set properties of potential
+uL  = 1e3 * phys.pc             ; unit of length
+uM  = phys.msun                 ; unit of mass
+uV  = 1d5                       ; unit of velocity
+
+; derived units
+constG   = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
+pcentre  = [50.,50.,50.] * 1d3 * pc / uL
+mextern  = 1d10 * msun / uM
+;
+;
+;
+ifile  = 0
+for ifile=0,nfiles-1 do begin
+;for ifile=0,3 do begin
+   inf    = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
+   time   = h5ra(inf, 'Header','Time')
+   p      = h5rd(inf,'PartType1/Coordinates')
+   v      = h5rd(inf,'PartType1/Velocities')
+   id     = h5rd(inf,'PartType1/ParticleIDs')
+   indx   = sort(id)
+;
+   id     = id[indx]
+   for ic=0,2 do begin
+      tmp = reform(p[ic,*]) & p[ic,*] = tmp[indx]
+      tmp = reform(v[ic,*]) & v[ic,*] = tmp[indx]
+   endfor
+; calculate energy
+   dd  = size(p,/dimen) & npart = dd[1]
+   ener = fltarr(npart)
+   dr   = fltarr(npart) & dv = dr
+   for ic=0,2 do dr[*] = dr[*] + (p[ic,*]-pcentre[ic])^2
+   for ic=0,2 do dv[*] = dv[*] + v[ic,*]^2
+   dr = sqrt(dr)
+;   print,'time = ',time,p[0,0],v[0,0],id[0]
+   ek   = 0.5 * dv
+   ep   = - constG * mextern / dr
+   ener = ek + ep
+   tout(ifile) = time
+   eout(*,ifile) = ener[0:nfollow-1]
+   ekin(*,ifile) = ek[0:nfollow-1]
+   epot(*,ifile) = ep[0:nfollow-1]
+endfor
+
+; calculate relative energy change
+de = 0.0 * eout
+for ifile=1, nfiles -1 do de[*,ifile] = (eout[*,ifile]-eout[*,0])/eout[*,0]
+
+
+end
+
+

examples/Makefile.am

@@ -62,21 +62,19 @@ swift_fixdt_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="e
 swift_fixdt_mpi_LDADD =  ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS)
 
 # Scripts to generate ICs
-EXTRA_DIST = UniformBox/makeIC.py \
+EXTRA_DIST = UniformBox/makeIC.py UniformBox/run.sh UniformBox/uniformBox.yml \
 	     UniformDMBox/makeIC.py \
 	     PerturbedBox/makeIC.py \
-	     SedovBlast/makeIC.py SedovBlast/makeIC_fcc.py SedovBlast/solution.py \
-	     SodShock/makeIC.py SodShock/solution.py SodShock/glass_001.hdf5 SodShock/glass_002.hdf5 SodShock/rhox.py \
-	     CosmoVolume/getIC.sh \
+	     SedovBlast/makeIC.py SedovBlast/makeIC_fcc.py SedovBlast/solution.py SedovBlast/run.sh SedovBlast/sedov.yml \
+	     SodShock/makeIC.py SodShock/solution.py SodShock/glass_001.hdf5 SodShock/glass_002.hdf5 SodShock/rhox.py SodShock/run.sh SodShock/sodShock.yml \
+	     CosmoVolume/getIC.sh CosmoVolume/run.sh CosmoVolume/cosmoVolume.yml \
 	     BigCosmoVolume/makeIC.py \
 	     BigPerturbedBox/makeIC_fcc.py \
              GreshoVortex/makeIC.py GreshoVortex/solution.py \
-             MultiTypes/makeIC.py
-
+             MultiTypes/makeIC.py \
+             parameter_example.yml
 
 # Scripts to plot task graphs
 EXTRA_DIST += plot_tasks_MPI.py plot_tasks.py \
 	      process_plot_tasks_MPI process_plot_tasks
 
-# Simple run scripts
-EXTRA_DIST += runs.sh

examples/SedovBlast/run.sh

@@ -7,4 +7,4 @@ then
     python makeIC_fcc.py
 fi
 
-../swift -s sedov.yml
+../swift -s -t 16 sedov.yml

examples/SedovBlast/sedov.yml

 
 # Define the system of units to use internally. 
-UnitSystem:
+InternalUnitSystem:
   UnitMass_in_cgs:     1   # Grams
   UnitLength_in_cgs:   1   # Centimeters
   UnitVelocity_in_cgs: 1   # Centimeters per second
@@ -9,7 +9,6 @@ UnitSystem:
 
 # 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).
@@ -22,13 +21,25 @@ TimeIntegration:
   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 governing the snapshots
+Snapshots:
+  basename:            sedov # Common part of the name of output files
+  time_first:          0.    # Time of the first output (in internal units)
+  delta_time:          0.1   # Time difference between consecutive outputs (in internal units)
+  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 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.
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
   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).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  max_volume_change:     2.       # Maximal allowed change of volume over one time-step
 
 # Parameters related to the initial conditions
 InitialConditions:

examples/SodShock/run.sh

@@ -7,4 +7,4 @@ then
     python makeIC.py
 fi
 
-../swift -s sodShock.yml
+../swift -s -t 16 sodShock.yml

examples/SodShock/sodShock.yml

 
 # Define the system of units to use internally. 
-UnitSystem:
+InternalUnitSystem:
   UnitMass_in_cgs:     1   # Grams
   UnitLength_in_cgs:   1   # Centimeters
   UnitVelocity_in_cgs: 1   # Centimeters per second
@@ -9,7 +9,6 @@ UnitSystem:
 
 # 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.
@@ -22,13 +21,25 @@ TimeIntegration:
   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 governing the snapshots
+Snapshots:
+  basename:            sod # Common part of the name of output files
+  time_first:          0.  # Time of the first output (in internal units)
+  delta_time:          0.1 # Time difference between consecutive outputs (in internal units)
+  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 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.
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
   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).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  max_volume_change:     2.       # Maximal allowed change of volume over one time-step
 
 # Parameters related to the initial conditions
 InitialConditions:

examples/UniformBox/run.sh

@@ -7,4 +7,4 @@ then
     python makeIC.py 100
 fi
 
-../swift -s uniformBox.yml
+../swift -s -t 16 uniformBox.yml

examples/UniformBox/uniformBox.yml

 
 # Define the system of units to use internally. 
-UnitSystem:
+InternalUnitSystem:
   UnitMass_in_cgs:     1   # Grams
   UnitLength_in_cgs:   1   # Centimeters
   UnitVelocity_in_cgs: 1   # Centimeters per second
@@ -9,7 +9,6 @@ UnitSystem:
 
 # 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).
@@ -22,14 +21,26 @@ TimeIntegration:
   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 governing the snapshots
+Snapshots:
+  basename:            uniformBox # Common part of the name of output files
+  time_first:          0.         # Time of the first output (in internal units)
+  delta_time:          0.01       # Time difference between consecutive outputs (in internal units)
+  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 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.
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
   max_ghost_iterations:  30       # Maximal number of iterations allowed to converge towards the smoothing length.
   max_smoothing_length:  0.1      # Maximal smoothing length allowed (in internal units).
-
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  max_volume_change:     2.       # Maximal allowed change of volume over one time-step
+  
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  ./uniformBox.hdf5     # The file to read

examples/main.c

@@ -3,6 +3,9 @@
  * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
  *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
  *               2015 Peter W. Draper (p.w.draper@durham.ac.uk)
+ *                    Angus Lepper (angus.lepper@ed.ac.uk)
+ *               2016 John A. Regan (john.a.regan@durham.ac.uk)
+ *                    Tom Theuns (tom.theuns@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
@@ -67,6 +70,9 @@ void print_help_message() {
          "Run with an external gravitational potential");
   printf("  %2s %8s %s\n", "-G", "", "Run with self-gravity");
   printf("  %2s %8s %s\n", "-s", "", "Run with SPH");
+  printf("  %2s %8s %s\n", "-t", "{int}",
+         "The number of threads to use on each MPI rank. Defaults to 1 if not "
+         "specified.");
   printf("  %2s %8s %s\n", "-v", "[12]",
          "Increase the level of verbosity 1: MPI-rank 0 writes ");
   printf("  %2s %8s %s\n", "", "", "2: All MPI-ranks write");
@@ -134,12 +140,13 @@ int main(int argc, char *argv[]) {
   int with_hydro = 0;
   int with_fp_exceptions = 0;
   int verbose = 0;
+  int nr_threads = 1;
   char paramFileName[200] = "";
   unsigned long long cpufreq = 0;
 
   /* Parse the parameters */
   int c;
-  while ((c = getopt(argc, argv, "cdef:gGhsv:y")) != -1) switch (c) {
+  while ((c = getopt(argc, argv, "cdef:gGhst:v:y:")) != -1) switch (c) {
       case 'c':
         with_cosmology = 1;
         break;
@@ -168,6 +175,14 @@ int main(int argc, char *argv[]) {
       case 's':
         with_hydro = 1;
         break;
+      case 't':
+        if (sscanf(optarg, "%d", &nr_threads) != 1) {
+          if (myrank == 0)
+            printf("Error parsing the number of threads (-t).\n");
+          if (myrank == 0) print_help_message();
+          return 1;
+        }
+        break;
       case 'v':
         if (sscanf(optarg, "%d", &verbose) != 1) {
           if (myrank == 0) printf("Error parsing verbosity level (-v).\n");
@@ -249,17 +264,6 @@ int main(int argc, char *argv[]) {
   MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
 #endif
 
-  /* Initialize unit system */
-  struct UnitSystem us;
-  units_init(&us, params);
-  if (myrank == 0) {
-    message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
-    message("Unit system: U_L = %e cm.", us.UnitLength_in_cgs);
-    message("Unit system: U_t = %e s.", us.UnitTime_in_cgs);
-    message("Unit system: U_I = %e A.", us.UnitCurrent_in_cgs);
-    message("Unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
-  }
-
 /* Prepare the domain decomposition scheme */
 #ifdef WITH_MPI
   struct partition initial_partition;
@@ -277,6 +281,29 @@ int main(int argc, char *argv[]) {
   }
 #endif
 
+  /* Initialize unit system and constants */
+  struct UnitSystem us;
+  struct phys_const prog_const;
+  units_init(&us, params, "InternalUnitSystem");
+  phys_const_init(&us, &prog_const);
+  if (myrank == 0 && verbose > 0) {
+    message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
+    message("Unit system: U_L = %e cm.", us.UnitLength_in_cgs);
+    message("Unit system: U_t = %e s.", us.UnitTime_in_cgs);
+    message("Unit system: U_I = %e A.", us.UnitCurrent_in_cgs);
+    message("Unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
+    phys_const_print(&prog_const);
+  }
+
+  /* Initialise the hydro properties */
+  struct hydro_props hydro_properties;
+  hydro_props_init(&hydro_properties, params);
+
+  /* Initialise the external potential properties */
+  struct external_potential potential;
+  if (with_external_gravity) potential_init(params, &us, &potential);
+  if (with_external_gravity && myrank == 0) potential_print(&potential);
+
   /* Read particles and space information from (GADGET) ICs */
   char ICfileName[200] = "";
   parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
@@ -314,6 +341,13 @@ int main(int argc, char *argv[]) {
     for (size_t k = 0; k < Ngas; ++k) parts[k].gpart = NULL;
     Ngpart = 0;
   }
+  if (!with_hydro) {
+    free(parts);
+    parts = NULL;
+    for (size_t k = 0; k < Ngpart; ++k)
+      if (gparts[k].id > 0) error("Linking problem");
+    Ngas = 0;
+  }
 
   /* Get the total number of particles across all nodes. */
   long long N_total[2] = {0, 0};
@@ -377,7 +411,8 @@ int main(int argc, char *argv[]) {
   /* Initialize the engine with the space and policies. */
   if (myrank == 0) clocks_gettime(&tic);
   struct engine e;
-  engine_init(&e, &s, params, nr_nodes, myrank, engine_policies, talking);
+  engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, engine_policies,
+              talking, &prog_const, &hydro_properties, &potential);
   if (myrank == 0) {
     clocks_gettime(&toc);
     message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
@@ -385,33 +420,13 @@ int main(int argc, char *argv[]) {
     fflush(stdout);
   }
 
+  /* Write the state of the system before starting time integration. */
+  if (!dry_run) engine_dump_snapshot(&e);
+
   /* Now that everything is ready, no need for the parameters any more */
   free(params);
   params = NULL;
 
-  int with_outputs = 1;
-  if (with_outputs && !dry_run) {
-    /* Write the state of the system before starting time integration. */
-    if (myrank == 0) clocks_gettime(&tic);
-#if defined(WITH_MPI)
-#if defined(HAVE_PARALLEL_HDF5)
-    write_output_parallel(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD,
-                          MPI_INFO_NULL);
-#else
-    write_output_serial(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD,
-                        MPI_INFO_NULL);
-#endif
-#else
-    write_output_single(&e, &us);
-#endif
-    if (myrank == 0 && verbose) {
-      clocks_gettime(&toc);
-      message("writing particle properties took %.3f %s.",
-              clocks_diff(&tic, &toc), clocks_getunit());
-      fflush(stdout);
-    }
-  }
-
 /* Init the runner history. */
 #ifdef HIST
   for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
@@ -420,10 +435,10 @@ int main(int argc, char *argv[]) {
   /* Get some info to the user. */
   if (myrank == 0) {
     message(
-        "Running on %lld gas particles and %lld DM particles until t=%.3e with "
-        "%i threads and %i queues (dt_min=%.3e, dt_max=%.3e)...",
-        N_total[0], N_total[1], e.timeEnd, e.nr_threads, e.sched.nr_queues,
-        e.dt_min, e.dt_max);
+        "Running on %lld gas particles and %lld DM particles from t=%.3e until "
+        "t=%.3e with %d threads and %d queues (dt_min=%.3e, dt_max=%.3e)...",
+        N_total[0], N_total[1], e.timeBegin, e.timeEnd, e.nr_threads,
+        e.sched.nr_queues, e.dt_min, e.dt_max);
     fflush(stdout);
   }
 
@@ -468,28 +483,6 @@ int main(int argc, char *argv[]) {
     /* Take a step. */
     engine_step(&e);
 
-    if (with_outputs && j % 100 == 0) {
-
-      if (myrank == 0) clocks_gettime(&tic);
-#if defined(WITH_MPI)
-#if defined(HAVE_PARALLEL_HDF5)
-      write_output_parallel(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD,
-                            MPI_INFO_NULL);
-#else
-      write_output_serial(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD,
-                          MPI_INFO_NULL);
-#endif
-#else
-      write_output_single(&e, &us);
-#endif
-      if (myrank == 0 && verbose) {
-        clocks_gettime(&toc);
-        message("writing particle properties took %.3f %s.",
-                clocks_diff(&tic, &toc), clocks_getunit());
-        fflush(stdout);
-      }
-    }
-
     /* Dump the task data using the given frequency. */
     if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
 #ifdef WITH_MPI
@@ -571,28 +564,8 @@ int main(int argc, char *argv[]) {
            (double)runner_hist_bins[k]);
 #endif
 
-  if (with_outputs) {
-
-    if (myrank == 0) clocks_gettime(&tic);
-/* Write final output. */
-#if defined(WITH_MPI)
-#if defined(HAVE_PARALLEL_HDF5)
-    write_output_parallel(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD,
-                          MPI_INFO_NULL);
-#else
-    write_output_serial(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD,
-                        MPI_INFO_NULL);
-#endif
-#else
-    write_output_single(&e, &us);
-#endif
-    if (myrank == 0 && verbose) {
-      clocks_gettime(&toc);
-      message("writing particle properties took %.3f %s.",
-              clocks_diff(&tic, &toc), clocks_getunit());
-      fflush(stdout);
-    }
-  }
+  /* Write final output. */
+  engine_dump_snapshot(&e);
 
 #ifdef WITH_MPI
   if ((res = MPI_Finalize()) != MPI_SUCCESS)

examples/parameter_example.yml

 
 # Define the system of units to use internally. 
-UnitSystem:
+InternalUnitSystem:
   UnitMass_in_cgs:     1   # Grams
   UnitLength_in_cgs:   1   # Centimeters
   UnitVelocity_in_cgs: 1   # Centimeters per second
@@ -9,7 +9,6 @@ UnitSystem:
 
 # Parameters for the task scheduling
 Scheduler:
-  nr_threads:       2        # 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:    8000000  # Maximal number of interactions per sub-task  (this is the default value).
@@ -22,13 +21,25 @@ TimeIntegration:
   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 governing the snapshots
+Snapshots:
+  basename:   output      # Common part of the name of output files
+  time_first: 0.          # Time of the first output (in internal units)
+  delta_time: 0.01        # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1  # Unit system for the outputs (Grams)
+  UnitLength_in_cgs:   1  # Unit system for the outputs (Centimeters)
+  UnitVelocity_in_cgs: 1  # Unit system for the outputs (Centimeters per second)
+  UnitCurrent_in_cgs:  1  # Unit system for the outputs (Amperes)
+  UnitTemp_in_cgs:     1  # Unit system for the outputs (Kelvin)
+
 # 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.
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
   max_ghost_iterations:  30       # Maximal number of iterations allowed to converge towards the smoothing length.
-  max_smoothing_length:  3.       # Maximal smoothing length allowed (in internal units).
+  max_smoothing_length:  0.1      # Maximal smoothing length allowed (in internal units).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  max_volume_change:     2.       # Maximal allowed change of volume over one time-step
 
 # Parameters related to the initial conditions
 InitialConditions:
@@ -46,3 +57,11 @@ DomainDecomposition:
   initial_grid_z:    10
   repartition_type:   b     # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
  
+# Parameters related to external potentials
+  
+# Point mass external potential
+PointMass:
+  position_x:      50.     # location of external point mass in internal units
+  position_y:      50.
+  position_z:      50.
+  mass:            1e10     # mass of external point mass in internal units

examples/plot_tasks.py

@@ -58,9 +58,10 @@ PLOT_PARAMS = {"axes.labelsize": 10,
 pl.rcParams.update(PLOT_PARAMS)
 
 #  Tasks and subtypes. Indexed as in tasks.h.
-TASKTYPES = ["none", "sort", "self", "pair", "sub", "init", "ghost", "drift", "kick",
-             "send", "recv", "grav_pp", "grav_mm", "grav_up", "grav_down",
-             "part_sort", "gpart_sort", "split_cell", "rewait", "count"]
+TASKTYPES = ["none", "sort", "self", "pair", "sub", "init", "ghost", "drift",
+             "kick", "send", "recv", "grav_pp", "grav_mm", "grav_up",
+             "grav_down", "grav_external", "part_sort", "gpart_sort",
+             "split_cell", "rewait", "count"]
 
 TASKCOLOURS = {"none": "black",
                "sort": "lightblue",
@@ -77,6 +78,7 @@ TASKCOLOURS = {"none": "black",
                "grav_mm": "mediumturquoise",
                "grav_up": "mediumvioletred",
                "grav_down": "mediumnightblue",
+               "grav_external": "darkred",
                "part_sort": "steelblue",
                "gpart_sort": "teal" ,
                "split_cell": "seagreen",

examples/plot_tasks_MPI.py

@@ -64,9 +64,10 @@ PLOT_PARAMS = {"axes.labelsize": 10,
 pl.rcParams.update(PLOT_PARAMS)
 
 #  Tasks and subtypes. Indexed as in tasks.h.
-TASKTYPES = ["none", "sort", "self", "pair", "sub", "init", "ghost", "drift", "kick",
-             "send", "recv", "grav_pp", "grav_mm", "grav_up", "grav_down",
-             "part_sort", "gpart_sort", "split_cell", "rewait", "count"]
+TASKTYPES = ["none", "sort", "self", "pair", "sub", "init", "ghost", "drift",
+             "kick", "send", "recv", "grav_pp", "grav_mm", "grav_up",
+             "grav_down", "grav_external", "part_sort", "gpart_sort",
+             "split_cell", "rewait", "count"]
 
 TASKCOLOURS = {"none": "black",
                "sort": "lightblue",
@@ -83,8 +84,9 @@ TASKCOLOURS = {"none": "black",
                "grav_mm": "mediumturquoise",
                "grav_up": "mediumvioletred",
                "grav_down": "mediumnightblue",
+               "grav_external": "darkred",
                "part_sort": "steelblue",
-               "gpart_sort": "teal",
+               "gpart_sort": "teal" ,
                "split_cell": "seagreen",
                "rewait": "olive",
                "count": "powerblue"}