Merge branch 'master' into rebuild_criteria

Conflicts:
	src/cell.c
	src/cell.h
	src/runner.c
	src/runner_doiact.h

.gitignore

@@ -22,9 +22,6 @@ doc/Doxyfile
 
 examples/swift
 examples/swift_mpi
-examples/*.xmf
-examples/used_parameters.yml
-examples/energy.txt
 examples/*/*.xmf
 examples/*/*.hdf5
 examples/*/*.png

INSTALL.swift

@@ -143,6 +143,8 @@ before you can build it.
  - DOXYGEN: the doxygen library is required to create the SWIFT API
             documentation.
 
+ - python:  Examples and solution script use python and rely on the
+   	    numpy library version 1.8.2 or higher.
 
 
                              SWIFT Coding style

README

@@ -15,24 +15,26 @@ Usage: swift [OPTION]... PARAMFILE
        swift_mpi [OPTION]... PARAMFILE
 
 Valid options are:
-  -a          Pin runners using processor affinity
-  -c          Run with cosmological time integration
-  -C          Run with cooling
+  -a          Pin runners using processor affinity.
+  -c          Run with cosmological time integration.
+  -C          Run with cooling.
   -d          Dry run. Read the parameter file, allocate memory but does not read 
               the particles from ICs and exit before the start of time integration.
               Allows user to check validy of parameter and IC files as well as memory limits.
-  -D          Always drift all particles even the ones far from active particles.
-  -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
+  -D          Always drift all particles even the ones far from active particles. This emulates
+  	      Gadget-[23] and GIZMO's default behaviours.
+  -e          Enable floating-point exceptions (debugging mode).
+  -f    {int} Overwrite the CPU frequency (Hz) to be used for time measurements.
+  -g          Run with an external gravitational potential.
+  -G          Run with self-gravity.
   -n    {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop. 
-  -s          Run with SPH
+  -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.
   -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
+  -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

@@ -16,7 +16,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 # Init the project.
-AC_INIT([SWIFT],[0.4.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
+AC_INIT([SWIFT],[0.5.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
 swift_config_flags="$*"
 
 #  Need to define this, instead of using fifth argument of AC_INIT, until 2.64.
@@ -105,8 +105,9 @@ if test "$enable_mpi" = "yes"; then
        AC_PATH_PROG([MPIRUN],[mpirun],[notfound])
     fi
     if test "$MPIRUN" = "notfound"; then
+       # This may not be fatal (some systems do not allow mpirun on
+       # development nodes)., so push on.
        AC_MSG_WARN([Cannot find mpirun command on PATH, thread support may not be correct])
-       enable_mpi="no"
     else
        # Special options we know about.
        # Intel: -mt_mpi
@@ -200,6 +201,20 @@ if test "$enable_debugging_checks" = "yes"; then
    AC_DEFINE([SWIFT_DEBUG_CHECKS],1,[Enable expensive debugging])
 fi
 
+# Check if gravity force checks are on for some particles.
+AC_ARG_ENABLE([gravity-force-checks],
+   [AS_HELP_STRING([--enable-gravity-force-checks],
+     [Activate expensive brute-force gravity checks for a fraction 1/N of all particles @<:@N@:>@]
+   )],
+   [gravity_force_checks="$enableval"],
+   [gravity_force_checks="no"]
+)
+if test "$gravity_force_checks" == "yes"; then
+   AC_MSG_ERROR(Need to specify the fraction of particles to check when using --enable-gravity-force-checks!)
+elif test "$gravity_force_checks" != "no"; then
+   AC_DEFINE_UNQUOTED([SWIFT_GRAVITY_FORCE_CHECKS], [$enableval] ,[Enable gravity brute-force checks])
+fi
+
 # Define HAVE_POSIX_MEMALIGN if it works.
 AX_FUNC_POSIX_MEMALIGN
 
@@ -637,13 +652,13 @@ AC_ARG_WITH([hydro-dimension],
 )
 case "$with_dimension" in
    1)
-      AC_DEFINE([HYDRO_DIMENSION_1D], [1], [1D analysis])
+      AC_DEFINE([HYDRO_DIMENSION_1D], [1], [1D solver])
    ;; 
    2)
-      AC_DEFINE([HYDRO_DIMENSION_2D], [2], [2D analysis])
+      AC_DEFINE([HYDRO_DIMENSION_2D], [2], [2D solver])
    ;; 
    3)
-      AC_DEFINE([HYDRO_DIMENSION_3D], [3], [3D analysis])
+      AC_DEFINE([HYDRO_DIMENSION_3D], [3], [3D solver])
    ;; 
    *)
       AC_MSG_ERROR([Dimensionality must be 1, 2 or 3])
@@ -751,7 +766,7 @@ esac
 #  External potential
 AC_ARG_WITH([ext-potential],
    [AS_HELP_STRING([--with-ext-potential=<pot>],
-      [external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch default: none@:>@]
+      [external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch, sine-wave default: none@:>@]
    )],
    [with_potential="$withval"],
    [with_potential="none"]
@@ -769,6 +784,9 @@ case "$with_potential" in
    disc-patch)
       AC_DEFINE([EXTERNAL_POTENTIAL_DISC_PATCH], [1], [Disc-patch external potential])
    ;; 
+   sine-wave)
+      AC_DEFINE([EXTERNAL_POTENTIAL_SINE_WAVE], [1], [Sine wave external potential in 1D])
+   ;; 
    *)
       AC_MSG_ERROR([Unknown external potential: $with_potential])
    ;;
@@ -821,8 +839,10 @@ AC_MSG_RESULT([
    Riemann solver     : $with_riemann
    Cooling function   : $with_cooling
    External potential : $with_potential
+
    Task debugging     : $enable_task_debugging
    Debugging checks   : $enable_debugging_checks
+   Gravity checks     : $gravity_force_checks
 ])
 
 # Make sure the latest git revision string gets included

doc/Doxyfile.in

@@ -762,6 +762,7 @@ WARN_LOGFILE           =
 INPUT                  =  @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
 INPUT		       += @top_srcdir@/src/hydro/Minimal
 INPUT		       += @top_srcdir@/src/gravity/Default
+INPUT		       += @top_srcdir@/src/stars/Default
 INPUT		       += @top_srcdir@/src/riemann
 INPUT		       += @top_srcdir@/src/potential/point_mass
 INPUT		       += @top_srcdir@/src/cooling/const_du

examples/BigCosmoVolume/makeIC.py

@@ -107,7 +107,7 @@ if n_copy > 1:
     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, j * boxSize, k * boxSize ]), axis=0)
+                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)

examples/CoolingHalo/makeIC.py

@@ -227,7 +227,7 @@ ds[()] = u
 u = np.zeros(1)
 
 # Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False, dtype='L')
+ids = 1 + np.linspace(0, N, N, endpoint=False)
 ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
 ds[()] = ids
 

examples/CoolingHaloWithSpin/makeIC.py

@@ -233,7 +233,7 @@ ds[()] = u
 u = np.zeros(1)
 
 # Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False, dtype='L')
+ids = 1 + np.linspace(0, N, N, endpoint=False)
 ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
 ds[()] = ids
 

examples/DiscPatch/GravityOnly/makeIC.py

@@ -150,7 +150,7 @@ ds[()] = m
 m = numpy.zeros(1)
 
 
-ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False, dtype='L')
+ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False)
 ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
 ds[()] = ids
 

examples/DiscPatch/HydroStatic/makeIC.py

@@ -205,7 +205,7 @@ if (entropy_flag == 1):
 else:
     ds[()] = u    
 
-ids = 1 + numpy.linspace(0, numGas, numGas, endpoint=False, dtype='L')
+ids = 1 + numpy.linspace(0, numGas, numGas, endpoint=False)
 ds = grp0.create_dataset('ParticleIDs', (numGas, ), 'L')
 ds[()] = ids
 

examples/ExternalPointMass/energy_plot.py

@@ -91,6 +91,9 @@ for i in range(402):
     E_tot_snap[i] = E_kin_snap[i] + E_pot_snap[i]
     Lz_snap[i] = np.sum(Lz)
 
+print "Starting energy:", E_kin_stats[0], E_pot_stats[0], E_tot_stats[0]
+print "Ending   energy:", E_kin_stats[-1], E_pot_stats[-1], E_tot_stats[-1]
+    
 # Plot energy evolution
 figure()
 plot(time_stats, E_kin_stats, "r-", lw=0.5, label="Kinetic energy")

examples/GreshoVortex_2D/plotSolution.py

@@ -31,6 +31,7 @@ P0 = 0.           # Constant additional pressure (should have no impact on the d
 import matplotlib
 matplotlib.use("Agg")
 from pylab import *
+from scipy import stats
 import h5py
 
 # Plot parameters
@@ -104,6 +105,26 @@ u = sim["/PartType0/InternalEnergy"][:]
 S = sim["/PartType0/Entropy"][:]
 P = sim["/PartType0/Pressure"][:]
 
+# Bin te data
+r_bin_edge = np.arange(0., 1., 0.02)
+r_bin = 0.5*(r_bin_edge[1:] + r_bin_edge[:-1])
+rho_bin,_,_ = stats.binned_statistic(r, rho, statistic='mean', bins=r_bin_edge)
+v_bin,_,_ = stats.binned_statistic(r, v_phi, statistic='mean', bins=r_bin_edge)
+P_bin,_,_ = stats.binned_statistic(r, P, statistic='mean', bins=r_bin_edge)
+S_bin,_,_ = stats.binned_statistic(r, S, statistic='mean', bins=r_bin_edge)
+u_bin,_,_ = stats.binned_statistic(r, u, statistic='mean', bins=r_bin_edge)
+rho2_bin,_,_ = stats.binned_statistic(r, rho**2, statistic='mean', bins=r_bin_edge)
+v2_bin,_,_ = stats.binned_statistic(r, v_phi**2, statistic='mean', bins=r_bin_edge)
+P2_bin,_,_ = stats.binned_statistic(r, P**2, statistic='mean', bins=r_bin_edge)
+S2_bin,_,_ = stats.binned_statistic(r, S**2, statistic='mean', bins=r_bin_edge)
+u2_bin,_,_ = stats.binned_statistic(r, u**2, statistic='mean', bins=r_bin_edge)
+rho_sigma_bin = np.sqrt(rho2_bin - rho_bin**2)
+v_sigma_bin = np.sqrt(v2_bin - v_bin**2)
+P_sigma_bin = np.sqrt(P2_bin - P_bin**2)
+S_sigma_bin = np.sqrt(S2_bin - S_bin**2)
+u_sigma_bin = np.sqrt(u2_bin - u_bin**2)
+
+
 # Plot the interesting quantities
 figure()
 
@@ -113,6 +134,7 @@ subplot(231)
 
 plot(r, v_phi, '.', color='r', ms=0.5)
 plot(solution_r, solution_v_phi, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, v_bin, yerr=v_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
 plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
@@ -125,6 +147,7 @@ subplot(232)
 
 plot(r, rho, '.', color='r', ms=0.5)
 plot(solution_r, solution_rho, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, rho_bin, yerr=rho_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
 plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
@@ -138,6 +161,7 @@ subplot(233)
 
 plot(r, P, '.', color='r', ms=0.5)
 plot(solution_r, solution_P, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, P_bin, yerr=P_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
 plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
@@ -150,6 +174,7 @@ subplot(234)
 
 plot(r, u, '.', color='r', ms=0.5)
 plot(solution_r, solution_u, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, u_bin, yerr=u_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
 plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)
@@ -163,6 +188,7 @@ subplot(235)
 
 plot(r, S, '.', color='r', ms=0.5)
 plot(solution_r, solution_s, '--', color='k', alpha=0.8, lw=1.2)
+errorbar(r_bin, S_bin, yerr=S_sigma_bin, fmt='.', ms=8.0, color='b', lw=1.2)
 plot([0.2, 0.2], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
 plot([0.4, 0.4], [-100, 100], ':', color='k', alpha=0.4, lw=1.2)
 xlabel("${\\rm{Radius}}~r$", labelpad=0)

examples/HydrostaticHalo/density_profile.py

@@ -28,7 +28,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
 #lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
 #X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])

examples/HydrostaticHalo/internal_energy_profile.py

@@ -46,7 +46,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
 #lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
 #X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])

examples/HydrostaticHalo/makeIC.py

@@ -227,7 +227,7 @@ ds[()] = u
 u = np.zeros(1)
 
 # Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False, dtype='L')
+ids = 1 + np.linspace(0, N, N, endpoint=False)
 ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
 ds[()] = ids
 

examples/HydrostaticHalo/velocity_profile.py

@@ -46,7 +46,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
 unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
 unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
 unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
 v_c_cgs = v_c * unit_velocity_cgs
 header = f["Header"]
 N = header.attrs["NumPart_Total"][0]

examples/IsothermalPotential/makeIC.py

@@ -138,7 +138,7 @@ ds = grp1.create_dataset('Masses', (numPart,), 'f')
 ds[()] = m
 m = numpy.zeros(1)
 
-ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False, dtype='L')
+ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False)
 ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
 ds[()] = ids
 

examples/IsothermalPotential/run.sh

@@ -4,7 +4,7 @@
 if [ ! -e Isothermal.hdf5 ]
 then
     echo "Generating initial conditions for the isothermal potential box example..."
-    python makeIC.py 1000 1
+    python makeIC.py 1000 0
 fi
 
 rm -rf Isothermal_*.hdf5

examples/MultiTypes/makeIC.py

@@ -36,110 +36,154 @@ eta = 1.2349             # 48 ngbs with cubic spline kernel
 rhoDM = 1.
 Ldm = int(sys.argv[2])  # Number of particles along one axis
 
-fileName = "multiTypes.hdf5"
+massStars = 0.1
+Lstars = int(sys.argv[3])  # Number of particles along one axis
+
+fileBaseName = "multiTypes"
+num_files = int(sys.argv[4])
 
 #---------------------------------------------------
-numGas = Lgas**3
-massGas = boxSize**3 * rhoGas / numGas
+numGas_tot = Lgas**3
+massGas = boxSize**3 * rhoGas / numGas_tot
 internalEnergy = P / ((gamma - 1.)*rhoGas)
 
-numDM = Ldm**3
-massDM = boxSize**3 * rhoDM / numDM
+numDM_tot = Ldm**3
+massDM = boxSize**3 * rhoDM / numDM_tot
+
+numStars_tot = Lstars**3
+massStars = massDM * massStars
+
 
 #--------------------------------------------------
 
-#File
-file = h5py.File(fileName, 'w')
-
-# Header
-grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] =  [numGas, numDM, 0, 0, 0, 0]
-grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [numGas, numDM, 0, 0, 0, 0]
-grp.attrs["Time"] = 0.0
-grp.attrs["NumFilesPerSnapshot"] = 1
-grp.attrs["MassTable"] = [0.0, massDM, 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"] = 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")
-
-v  = zeros((numGas, 3))
-ds = grp.create_dataset('Velocities', (numGas, 3), 'f')
-ds[()] = v
-v = zeros(1)
-
-m = full((numGas, 1), massGas)
-ds = grp.create_dataset('Masses', (numGas,1), 'f')
-ds[()] = m
-m = zeros(1)
-
-h = full((numGas, 1), eta * boxSize / Lgas)
-ds = grp.create_dataset('SmoothingLength', (numGas,1), 'f')
-ds[()] = h
-h = zeros(1)
-
-u = full((numGas, 1), internalEnergy)
-ds = grp.create_dataset('InternalEnergy', (numGas,1), 'f')
-ds[()] = u
-u = zeros(1)
-
-ids = linspace(0, numGas, numGas, endpoint=False).reshape((numGas,1))
-ds = grp.create_dataset('ParticleIDs', (numGas, 1), 'L')
-ds[()] = ids + 1
-x      = ids % Lgas;
-y      = ((ids - x) / Lgas) % Lgas;
-z      = (ids - x - Lgas * y) / Lgas**2;
-coords = zeros((numGas, 3))
-coords[:,0] = z[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
-coords[:,1] = y[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
-coords[:,2] = x[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
-ds = grp.create_dataset('Coordinates', (numGas, 3), 'd')
-ds[()] = coords
-
-
-
-
-
-# DM Particle group
-grp = file.create_group("/PartType1")
-
-v  = zeros((numDM, 3))
-ds = grp.create_dataset('Velocities', (numDM, 3), 'f')
-ds[()] = v
-v = zeros(1)
-
-m = full((numDM, 1), massDM)
-ds = grp.create_dataset('Masses', (numDM,1), 'f')
-ds[()] = m
-m = zeros(1)
-
-ids = linspace(0, numDM, numDM, endpoint=False).reshape((numDM,1))
-ds = grp.create_dataset('ParticleIDs', (numDM, 1), 'L')
-ds[()] = ids + Lgas**3 + 1
-x      = ids % Ldm;
-y      = ((ids - x) / Ldm) % Ldm;
-z      = (ids - x - Ldm * y) / Ldm**2;
-coords = zeros((numDM, 3))
-coords[:,0] = z[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
-coords[:,1] = y[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
-coords[:,2] = x[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
-ds = grp.create_dataset('Coordinates', (numDM, 3), 'd')
-ds[()] = coords
-
-file.close()
+offsetGas = 0
+offsetDM = 0
+offsetStars = 0
+
+for n in range(num_files):
+
+    # File name
+    if num_files == 1:
+        fileName = fileBaseName + ".hdf5"
+    else:
+        fileName = fileBaseName + ".%d.hdf5"%n
+        
+    # File
+    file = h5py.File(fileName, 'w')
+
+    # Number of particles
+    numGas = numGas_tot / num_files
+    numDM = numDM_tot / num_files
+    numStars = numStars_tot / num_files
+
+    if n == num_files - 1:
+        numGas += numGas_tot % num_files
+        numDM += numDM_tot % num_files
+        numStars += numStars_tot % num_files
+
+    
+    # Header
+    grp = file.create_group("/Header")
+    grp.attrs["BoxSize"] = boxSize
+    grp.attrs["NumPart_Total"] =  [numGas_tot, numDM_tot, 0, 0, numStars_tot, 0]
+    grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+    grp.attrs["NumPart_ThisFile"] = [numGas, numDM, 0, 0, numStars, 0]
+    grp.attrs["Time"] = 0.0
+    grp.attrs["NumFilesPerSnapshot"] = num_files
+    grp.attrs["MassTable"] = [0.0, massDM, 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"] = 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")
+
+    v  = zeros((numGas, 3))
+    ds = grp.create_dataset('Velocities', (numGas, 3), 'f', data=v)
+    
+    m = full((numGas, 1), massGas)
+    ds = grp.create_dataset('Masses', (numGas,1), 'f', data=m)
+    
+    h = full((numGas, 1), eta * boxSize / Lgas)
+    ds = grp.create_dataset('SmoothingLength', (numGas,1), 'f', data=h)
+    
+    u = full((numGas, 1), internalEnergy)
+    ds = grp.create_dataset('InternalEnergy', (numGas,1), 'f', data=u)
+
+    ids = linspace(offsetGas, offsetGas+numGas, numGas, endpoint=False).reshape((numGas,1))
+    ds = grp.create_dataset('ParticleIDs', (numGas, 1), 'L', data=ids+1)
+    x      = ids % Lgas;
+    y      = ((ids - x) / Lgas) % Lgas;
+    z      = (ids - x - Lgas * y) / Lgas**2;
+    coords = zeros((numGas, 3))
+    coords[:,0] = z[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
+    coords[:,1] = y[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
+    coords[:,2] = x[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
+    ds = grp.create_dataset('Coordinates', (numGas, 3), 'd', data=coords)
+
+
+    
+    # DM Particle group
+    grp = file.create_group("/PartType1")
+
+    v  = zeros((numDM, 3))
+    ds = grp.create_dataset('Velocities', (numDM, 3), 'f', data=v)
+
+    m = full((numDM, 1), massDM)
+    ds = grp.create_dataset('Masses', (numDM,1), 'f', data=m)
+
+    ids = linspace(offsetDM, offsetDM+numDM, numDM, endpoint=False).reshape((numDM,1))
+    ds = grp.create_dataset('ParticleIDs', (numDM, 1), 'L', data=ids + numGas_tot + 1)
+    ds[()] = ids + Lgas**3 + 1
+    x      = ids % Ldm;
+    y      = ((ids - x) / Ldm) % Ldm;
+    z      = (ids - x - Ldm * y) / Ldm**2;
+    coords = zeros((numDM, 3))
+    coords[:,0] = z[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
+    coords[:,1] = y[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
+    coords[:,2] = x[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
+    ds = grp.create_dataset('Coordinates', (numDM, 3), 'd', data=coords)
+
+
+
+    # Star Particle group
+    grp = file.create_group("/PartType4")
+
+    v  = zeros((numStars, 3))
+    ds = grp.create_dataset('Velocities', (numStars, 3), 'f', data=v)
+
+    m = full((numStars, 1), massStars)
+    ds = grp.create_dataset('Masses', (numStars,1), 'f', data=m)
+
+    ids = linspace(0, numStars, numStars, endpoint=False).reshape((numStars,1))
+    ds = grp.create_dataset('ParticleIDs', (numStars, 1), 'L', data=ids + numGas_tot + numDM_tot + 1)
+    x      = ids % Ldm;
+    y      = ((ids - x) / Ldm) % Ldm;
+    z      = (ids - x - Ldm * y) / Ldm**2;
+    coords = zeros((numStars, 3))
+    coords[:,0] = z[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
+    coords[:,1] = y[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
+    coords[:,2] = x[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
+    ds = grp.create_dataset('Coordinates', (numStars, 3), 'd', data=coords)
+
+
+    
+    # Shift stuff
+    offsetGas += numGas
+    offsetDM += numDM
+    offsetStars += numStars
+    
+    file.close()
+

examples/MultiTypes/multiTypes.yml

@@ -32,6 +32,7 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  ./multiTypes.hdf5     # The file to read
+  replicate:  2                     # Replicate all particles twice along each axis
 
 # External potential parameters
 PointMassPotential: