Merge branch 'master' into gizmo_set_internal_energy

.gitignore

@@ -80,9 +80,12 @@ tests/testRiemannHLLC
 tests/testMatrixInversion
 
 theory/latex/swift.pdf
-theory/kernel/kernels.pdf
-theory/kernel/kernel_derivatives.pdf
-theory/kernel/kernel_definitions.pdf
+theory/SPH/Kernels/kernels.pdf
+theory/SPH/Kernels/kernel_derivatives.pdf
+theory/SPH/Kernels/kernel_definitions.pdf
+theory/SPH/Flavours/sph_flavours.pdf
+theory/SPH/EoS/eos.pdf
+theory/SPH/*.pdf
 theory/paper_pasc/pasc_paper.pdf
 
 m4/libtool.m4

configure.ac

@@ -155,6 +155,21 @@ LT_INIT
 AC_PROG_CC_C99
 AC_C_INLINE
 
+# If debugging try to show inlined functions.
+if test "x$enable_debug" = "xyes"; then
+   #  Show inlined functions.
+   if test "$ax_cv_c_compiler_vendor" = "gnu"; then
+      # Would like to use -gdwarf and let the compiler pick a good version
+      # but that doesn't always work.
+      AX_CHECK_COMPILE_FLAG([-gdwarf -fvar-tracking-assignments],
+        [inline_EXTRA_FLAGS="-gdwarf -fvar-tracking-assignments"],
+        [inline_EXTRA_FLAGS="-gdwarf-2 -fvar-tracking-assignments"])
+      CFLAGS="$CFLAGS $inline_EXTRA_FLAGS"
+   elif test "$ax_cv_c_compiler_vendor" = "intel"; then
+      CFLAGS="$CFLAGS -debug inline-debug-info"
+   fi
+fi
+
 # Define HAVE_POSIX_MEMALIGN if it works.
 AX_FUNC_POSIX_MEMALIGN
 
@@ -451,7 +466,7 @@ AC_LINK_IFELSE([AC_LANG_PROGRAM(
 AC_MSG_RESULT($rtc_ok)
 
 # Add warning flags by default, if these can be used. Option =error adds
-# -Werror to GCC, clang and Intel.  Note do this last as compiler tests may 
+# -Werror to GCC, clang and Intel.  Note do this last as compiler tests may
 # become errors, if that's an issue don't use CFLAGS for these, use an AC_SUBST().
 AC_ARG_ENABLE([compiler-warnings],
    [AS_HELP_STRING([--enable-compiler-warnings],
@@ -461,7 +476,7 @@ AC_ARG_ENABLE([compiler-warnings],
    [enable_warn="error"]
 )
 if test "$enable_warn" != "no"; then
-   
+
     # AX_CFLAGS_WARN_ALL does not give good warning flags for the Intel compiler
     # We will do this by hand instead and only default to the macro for unknown compilers
     case "$ax_cv_c_compiler_vendor" in
@@ -475,7 +490,7 @@ if test "$enable_warn" != "no"; then
 	     AX_CFLAGS_WARN_ALL
 	  ;;
     esac
-    
+
     # Add a "choke on warning" flag if it exists
     if test "$enable_warn" = "error"; then
        case "$ax_cv_c_compiler_vendor" in

examples/DiscPatch/HydroStatic/dynamic.pro

@@ -8,7 +8,8 @@ iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare f
 @physunits
 
 indir    = './'
-basefile = 'Disk-Patch-dynamic_'
+;basefile = 'Disc-Patch-dynamic_'
+basefile = 'Disc-Patch_'
 
 ; set properties of potential
 uL   = phys.pc                  ; unit of length
@@ -16,18 +17,27 @@ uM   = phys.msun                ; unit of mass
 uV   = 1d5                      ; unit of velocity
 
 ; properties of patch
-surface_density = 10.
+surface_density = 100.          ; surface density of all mass, which generates the gravitational potential
 scale_height    = 100.
-z_disk          = 200.;
+z_disk          = 200.          ;
+fgas            = 0.1           ; gas fraction
 gamma           = 5./3.
 
 ; derived units
 constG   = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
 pcentre  = [0.,0.,z_disk] * pc / uL
 utherm     = !pi * constG * surface_density * scale_height / (gamma-1.)
+temp       = (utherm*uV^2)*phys.m_h/phys.kb
 soundspeed = sqrt(gamma * (gamma-1.) * utherm)
 t_dyn      = sqrt(scale_height / (constG * surface_density))
-
+rho0       = fgas*(surface_density)/(2.*scale_height)
+print,' dynamical time = ',t_dyn,' = ',t_dyn*UL/uV/(1d6*phys.yr),' Myr'
+print,' thermal energy per unit mass = ',utherm
+print,' central density = ',rho0,' = ',rho0*uM/uL^3/m_h,' particles/cm^3'
+print,' central temperature = ',temp
+lambda = 2 * !pi * phys.G^1.5 * (scale_height*uL)^1.5 * (surface_density * uM/uL^2)^0.5 * phys.m_h^2 / (gamma-1) / fgas
+print,' lambda = ',lambda
+stop
 ;
 infile = indir + basefile + '*'
 spawn,'ls -1 '+infile,res

examples/DiscPatch/HydroStatic/makeIC.py

@@ -56,9 +56,10 @@ print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
 
 
 # parameters of potential
-surface_density = 10.
+surface_density = 100. # surface density of all mass, which generates the gravitational potential
 scale_height    = 100.
 gamma           = 5./3.
+fgas            = 0.1  # gas fraction
 
 # derived units
 const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
@@ -131,7 +132,7 @@ h      = glass_h[0:numGas]
 numGas = numpy.shape(pos)[0]
 
 # compute furthe properties of ICs
-column_density = surface_density * numpy.tanh(boxSize/2./scale_height)
+column_density = fgas * surface_density * numpy.tanh(boxSize/2./scale_height)
 enclosed_mass  = column_density * boxSize * boxSize
 pmass          = enclosed_mass / numGas
 meanrho        = enclosed_mass / boxSize**3

examples/EAGLE_25/README

@@ -13,4 +13,4 @@ The particle load of the main EAGLE simulation can be reproduced by
 running these ICs on 64 cores.
 
 MD5 checksum of the ICs: 
-ada2c728db2bd2d77a20c4eef52dfaf1  EAGLE_ICs_25.hdf5
+02cd1c353b86230af047b5d4ab22afcf  EAGLE_ICs_25.hdf5

examples/Feedback/feedback.pro

+base = 'Feedback'
+inf  = 'Feedback_005.hdf5'
+
+blast  = [5.650488e-01, 5.004371e-01, 5.010494e-01] ; location of blast
+pos    = h5rd(inf,'PartType0/Coordinates')
+vel    = h5rd(inf,'PartType0/Velocities')
+rho    = h5rd(inf,'PartType0/Density')
+utherm = h5rd(inf,'PartType0/InternalEnergy')
+
+; shift to centre
+for ic=0,2 do pos[ic,*] = pos[ic,*] - blast[ic]
+
+;; distance from centre
+dist = fltarr(n_elements(rho))
+for ic=0,2 do dist = dist + pos[ic,*]^2
+dist = sqrt(dist)
+
+; radial velocity
+vr = fltarr(n_elements(rho))
+for ic=0,2 do vr = vr + pos[ic,*]*vel[ic,*]
+vr = vr / dist
+
+; 
+end

examples/Feedback/feedback.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  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 governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   5e-2  # 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-4  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            Feedback # 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)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-3 # Time between statistics output
+
+# 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_smoothing_length:  0.1       # Maximal smoothing length allowed (in internal units).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./Feedback.hdf5          # The file to read
+
+# Parameters for feedback
+
+SN:
+	time:    0.001 # time the SN explodes (internal units)
+	energy:  1.0   # energy of the explosion (internal units)
+	x:       0.5   # x-position of explostion (internal units)
+	y:       0.5   # y-position of explostion (internal units)
+	z:       0.5   # z-position of explostion (internal units)

examples/Feedback/makeIC.py

+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #               2016 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
+from numpy import *
+
+# Generates a swift IC file containing a cartesian distribution of particles
+# at a constant density and pressure in a cubic box
+
+# Parameters
+periodic= 1           # 1 For periodic box
+boxSize = 1.
+L = int(sys.argv[1])  # Number of particles along one axis
+rho = 1.              # Density
+P = 1.e-6             # Pressure
+gamma = 5./3.         # Gas adiabatic index
+eta = 1.2349          # 48 ngbs with cubic spline kernel
+fileName = "Feedback.hdf5" 
+
+#---------------------------------------------------
+numPart = L**3
+mass = boxSize**3 * rho / numPart
+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"] =  [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
+
+#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")
+
+v  = zeros((numPart, 3))
+ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+v = zeros(1)
+
+m = full((numPart, 1), mass)
+ds = grp.create_dataset('Masses', (numPart,1), 'f')
+ds[()] = m
+m = zeros(1)
+
+h = full((numPart, 1), eta * boxSize / L)
+ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
+ds[()] = h
+h = zeros(1)
+
+u = full((numPart, 1), internalEnergy)
+ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
+ds[()] = u
+u = zeros(1)
+
+
+ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
+ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
+ds[()] = ids + 1
+x      = ids % L;
+y      = ((ids - x) / L) % L;
+z      = (ids - x - L * y) / L**2;
+coords = zeros((numPart, 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 = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = coords
+
+file.close()

examples/SedovBlast_3D/plotSolution.py

@@ -18,7 +18,7 @@
  # 
  ##############################################################################
 
-# Computes the analytical solution of the 2D Sedov blast wave.
+# Computes the analytical solution of the 3D Sedov blast wave.
 # The script works for a given initial box and dumped energy and computes the solution at a later time t.
 
 # Parameters

examples/main.c

@@ -76,6 +76,7 @@ void print_help_message() {
          "Overwrite the CPU frequency (Hz) to be used for time measurements");
   printf("  %2s %8s %s\n", "-g", "",
          "Run with an external gravitational potential");
+  printf("  %2s %8s %s\n", "-F", "", "Run with feedback ");
   printf("  %2s %8s %s\n", "-G", "", "Run with self-gravity");
   printf("  %2s %8s %s\n", "-n", "{int}",
          "Execute a fixed number of time steps. When unset use the time_end "
@@ -147,6 +148,7 @@ int main(int argc, char *argv[]) {
   int nsteps = -2;
   int with_cosmology = 0;
   int with_external_gravity = 0;
+  int with_sourceterms = 0;
   int with_cooling = 0;
   int with_self_gravity = 0;
   int with_hydro = 0;
@@ -159,7 +161,7 @@ int main(int argc, char *argv[]) {
 
   /* Parse the parameters */
   int c;
-  while ((c = getopt(argc, argv, "acCdDef:gGhn:st:v:y:")) != -1) switch (c) {
+  while ((c = getopt(argc, argv, "acCdDef:FgGhn:st:v:y:")) != -1) switch (c) {
       case 'a':
         with_aff = 1;
         break;
@@ -185,6 +187,9 @@ int main(int argc, char *argv[]) {
           return 1;
         }
         break;
+      case 'F':
+        with_sourceterms = 1;
+        break;
       case 'g':
         with_external_gravity = 1;
         break;
@@ -257,15 +262,24 @@ int main(int argc, char *argv[]) {
     message(
         "Executing a dry run. No i/o or time integration will be performed.");
 
-  /* Report CPU frequency. */
+  /* Report CPU frequency.*/
   cpufreq = clocks_get_cpufreq();
   if (myrank == 0) {
     message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
   }
 
+/* Report host name(s). */
+#ifdef WITH_MPI
+  if (myrank == 0 || verbose > 1) {
+    message("Rank %d running on: %s", myrank, hostname());
+  }
+#else
+  message("Running on: %s", hostname());
+#endif
+
   /* Do we choke on FP-exceptions ? */
   if (with_fp_exceptions) {
-    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW | FE_UNDERFLOW);
+    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
     if (myrank == 0) message("Floating point exceptions will be reported.");
   }
 
@@ -443,6 +457,11 @@ int main(int argc, char *argv[]) {
   if (with_cooling) cooling_init(params, &us, &prog_const, &cooling_func);
   if (with_cooling && myrank == 0) cooling_print(&cooling_func);
 
+  /* Initialise the feedback properties */
+  struct sourceterms sourceterms;
+  if (with_sourceterms) sourceterms_init(params, &us, &sourceterms);
+  if (with_sourceterms && myrank == 0) sourceterms_print(&sourceterms);
+
   /* Construct the engine policy */
   int engine_policies = ENGINE_POLICY | engine_policy_steal;
   if (with_drift_all) engine_policies |= engine_policy_drift_all;
@@ -451,13 +470,14 @@ int main(int argc, char *argv[]) {
   if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
   if (with_cosmology) engine_policies |= engine_policy_cosmology;
   if (with_cooling) engine_policies |= engine_policy_cooling;
+  if (with_sourceterms) engine_policies |= engine_policy_sourceterms;
 
   /* 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, nr_threads, with_aff,
               engine_policies, talking, &us, &prog_const, &hydro_properties,
-              &potential, &cooling_func);
+              &potential, &cooling_func, &sourceterms);
   if (myrank == 0) {
     clocks_gettime(&toc);
     message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),

examples/parameter_example.yml

@@ -12,6 +12,7 @@ Scheduler:
   cell_max_size:    8000000  # (Optional) Maximal number of interactions per task if we force the split (this is the default value).
   cell_sub_size:    64000000 # (Optional) Maximal number of interactions per sub-task  (this is the default value).
   cell_split_size:  400      # (Optional) Maximal number of particles per cell (this is the default value).
+  cell_max_count:   10000    # (Optional) Maximal number of particles per cell allowed before triggering a sanitizing (this is the default value).
 
 # Parameters governing the time integration
 TimeIntegration:

src/Makefile.am

@@ -43,7 +43,8 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
     common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h \
     physical_constants.h physical_constants_cgs.h potential.h version.h \
-    hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h
+    hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
+    sourceterms_struct.h
 
 
 # Common source files
@@ -52,7 +53,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     units.c common_io.c single_io.c multipole.c version.c map.c \
     kernel_hydro.c tools.c part.c partition.c clocks.c parser.c \
     physical_constants.c potential.c hydro_properties.c \
-    runner_doiact_fft.c threadpool.c cooling.c
+    runner_doiact_fft.c threadpool.c cooling.c sourceterms.c
 
 # Include files for distribution, not installation.
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
@@ -62,6 +63,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
 		 gravity.h gravity_io.h \
 		 gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
 		 gravity/Default/gravity_debug.h gravity/Default/gravity_part.h  \
+		 sourceterms.h \
 	 	 hydro.h hydro_io.h \
 		 hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
                  hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \

src/cell.c

@@ -679,6 +679,59 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset) {
     part_relink_parts(gparts, gcount, parts - parts_offset);
 }
 
+/**
+ * @brief Sanitizes the smoothing length values of cells by setting large
+ * outliers to more sensible values.
+ *
+ * We compute the mean and standard deviation of the smoothing lengths in
+ * logarithmic space and limit values to mean + 4 sigma.
+ *
+ * @param c The cell.
+ */
+void cell_sanitize(struct cell *c) {
+
+  const int count = c->count;
+  struct part *parts = c->parts;
+
+  /* First collect some statistics */
+  float h_mean = 0.f, h_mean2 = 0.f;
+  float h_min = FLT_MAX, h_max = 0.f;
+  for (int i = 0; i < count; ++i) {
+
+    const float h = logf(parts[i].h);
+    h_mean += h;
+    h_mean2 += h * h;
+    h_max = max(h_max, h);
+    h_min = min(h_min, h);
+  }
+  h_mean /= count;
+  h_mean2 /= count;
+  const float h_var = h_mean2 - h_mean * h_mean;
+  const float h_std = (h_var > 0.f) ? sqrtf(h_var) : 0.1f * h_mean;
+
+  /* Choose a cut */
+  const float h_limit = expf(h_mean + 4.f * h_std);
+
+  /* Be verbose this is not innocuous */
+  message("Cell properties: h_min= %f h_max= %f geometric mean= %f.",
+          expf(h_min), expf(h_max), expf(h_mean));
+
+  if (c->h_max > h_limit) {
+
+    message("Smoothing lengths will be limited to (mean + 4sigma)= %f.",
+            h_limit);
+
+    /* Apply the cut */
+    for (int i = 0; i < count; ++i) parts->h = min(parts[i].h, h_limit);
+
+    c->h_max = h_limit;
+
+  } else {
+
+    message("Smoothing lengths will not be limited.");
+  }
+}
+
 /**
  * @brief Converts hydro quantities to a valid state after the initial density
  * calculation

src/cell.h

@@ -168,6 +168,9 @@ struct cell {
   /* Task for cooling */
   struct task *cooling;
 
+  /* Task for source terms */
+  struct task *sourceterms;
+
   /* Number of tasks that are associated with this cell. */
   int nr_tasks;
 
@@ -218,6 +221,7 @@ struct cell {
 
 /* Function prototypes. */
 void cell_split(struct cell *c, ptrdiff_t parts_offset);
+void cell_sanitize(struct cell *c);
 int cell_locktree(struct cell *c);
 void cell_unlocktree(struct cell *c);
 int cell_glocktree(struct cell *c);

src/const.h

@@ -96,6 +96,10 @@
 //#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
 //#define EXTERNAL_POTENTIAL_DISC_PATCH
 
+/* Source terms */
+#define SOURCETERMS_NONE
+//#define SOURCETERMS_SN_FEEDBACK
+
 /* Cooling properties */
 #define COOLING_NONE
 //#define COOLING_CONST_DU

src/engine.c

@@ -68,7 +68,7 @@
 #include "units.h"
 #include "version.h"
 
-const char *engine_policy_names[15] = {"none",
+const char *engine_policy_names[16] = {"none",
                                        "rand",
                                        "steal",
                                        "keep",
@@ -82,7 +82,8 @@ const char *engine_policy_names[15] = {"none",
                                        "external_gravity",
                                        "cosmology_integration",
                                        "drift_all",
-                                       "cooling"};
+                                       "cooling",
+                                       "sourceterms"};
 
 /** The rank of the engine as a global variable (for messages). */
 int engine_rank;
@@ -187,6 +188,8 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
   const int is_fixdt = (e->policy & engine_policy_fixdt) == engine_policy_fixdt;
   const int is_with_cooling =
       (e->policy & engine_policy_cooling) == engine_policy_cooling;
+  const int is_with_sourceterms =
+      (e->policy & engine_policy_sourceterms) == engine_policy_sourceterms;
 
   /* Is this the super-cell? */
   if (super == NULL && (c->density != NULL || (c->count > 0 && !c->split))) {
@@ -216,16 +219,18 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
       /* Generate the ghost task. */
       c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0, 0,
                                    c, NULL, 0);
-
 #ifdef EXTRA_HYDRO_LOOP
       /* Generate the extra ghost task. */
       c->extra_ghost = scheduler_addtask(s, task_type_extra_ghost,
                                          task_subtype_none, 0, 0, c, NULL, 0);
 #endif
-
       if (is_with_cooling)
         c->cooling = scheduler_addtask(s, task_type_cooling, task_subtype_none,
                                        0, 0, c, NULL, 0);
+      /* add source terms */
+      if (is_with_sourceterms)
+        c->sourceterms = scheduler_addtask(s, task_type_sourceterms,
+                                           task_subtype_none, 0, 0, c, NULL, 0);
     }
   }
 
@@ -1791,12 +1796,18 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
       scheduler_addunlock(sched, t->ci->init, t);
       scheduler_addunlock(sched, t, t->ci->kick);
     }
-
-    /* Cooling tasks should depend on kick and does not unlock anything since
-     it is the last task*/
+    /* Cooling tasks should depend on kick and unlock sourceterms */
     else if (t->type == task_type_cooling) {
       scheduler_addunlock(sched, t->ci->kick, t);
     }
+    /* source terms depend on cooling if performed, else on kick. It is the last
+       task */
+    else if (t->type == task_type_sourceterms) {
+      if (e->policy == engine_policy_cooling)
+        scheduler_addunlock(sched, t->ci->cooling, t);
+      else
+        scheduler_addunlock(sched, t->ci->kick, t);
+    }
   }
 }
 
@@ -2243,6 +2254,8 @@ void engine_rebuild(struct engine *e) {
   /* Re-build the space. */
   space_rebuild(e->s, 0.0, e->verbose);
 
+  if (e->ti_current == 0) space_sanitize(e->s);
+
 /* If in parallel, exchange the cell structure. */
 #ifdef WITH_MPI
   engine_exchange_cells(e);
@@ -2818,6 +2831,11 @@ void engine_step(struct engine *e) {
     mask |= 1 << task_type_cooling;
   }