Merge branch 'master' into initialisation_fixes

CONTRIBUTING.md

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

examples/plot_tasks.py

@@ -56,8 +56,8 @@ pl.rcParams.update(PLOT_PARAMS)
 
 #  Tasks and subtypes. Indexed as in tasks.h.
 TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair", "init", "ghost",
-             "drift", "kick", "kick_fixdt", "send", "recv", "grav_pp", "grav_mm",
-             "grav_up", "grav_down", "grav_external", "part_sort", "gpart_sort",
+             "drift", "kick", "kick_fixdt", "send", "recv", "grav_gather_m", "grav_fft",
+             "grav_mm", "grav_up", "grav_external", "part_sort", "gpart_sort",
              "split_cell", "rewait", "count"]
 
 TASKCOLOURS = {"none": "black",
@@ -73,10 +73,10 @@ TASKCOLOURS = {"none": "black",
                "kick_fixdt": "green",
                "send": "yellow",
                "recv": "magenta",
-               "grav_pp": "mediumorchid",
+               "grav_gather_m": "mediumorchid",
+               "grav_fft": "mediumnightblue",
                "grav_mm": "mediumturquoise",
                "grav_up": "mediumvioletred",
-               "grav_down": "mediumnightblue",
                "grav_external": "darkred",
                "part_sort": "steelblue",
                "gpart_sort": "teal" ,
@@ -84,12 +84,13 @@ TASKCOLOURS = {"none": "black",
                "rewait": "olive",
                "count": "powerblue"}
 
-SUBTYPES = ["none", "density", "force", "grav", "count"]
+SUBTYPES = ["none", "density", "force", "grav", "tend", "count"]
 
 SUBCOLOURS = {"none": "black",
               "density": "red",
               "force": "blue",
               "grav": "indigo",
+              "tend": "grey"
               "count": "purple"}
 
 #  Show docs if help is requested.

examples/plot_tasks_MPI.py

@@ -62,8 +62,8 @@ pl.rcParams.update(PLOT_PARAMS)
 
 #  Tasks and subtypes. Indexed as in tasks.h.
 TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair", "init", "ghost",
-             "drift", "kick", "kick_fixdt", "send", "recv", "grav_pp", "grav_mm",
-             "grav_up", "grav_down", "grav_external", "part_sort", "gpart_sort",
+             "drift", "kick", "kick_fixdt", "send", "recv", "grav_gather_m", "grav_fft",
+             "grav_mm", "grav_up", "grav_external", "part_sort", "gpart_sort",
              "split_cell", "rewait", "count"]
 
 TASKCOLOURS = {"none": "black",
@@ -79,10 +79,10 @@ TASKCOLOURS = {"none": "black",
                "kick_fixdt": "green",
                "send": "yellow",
                "recv": "magenta",
-               "grav_pp": "mediumorchid",
+               "grav_gather_m": "mediumorchid",
+               "grav_fft": "mediumnightblue",
                "grav_mm": "mediumturquoise",
                "grav_up": "mediumvioletred",
-               "grav_down": "mediumnightblue",
                "grav_external": "darkred",
                "part_sort": "steelblue",
                "gpart_sort": "teal" ,
@@ -90,12 +90,13 @@ TASKCOLOURS = {"none": "black",
                "rewait": "olive",
                "count": "powerblue"}
 
-SUBTYPES = ["none", "density", "force", "grav", "count"]
+SUBTYPES = ["none", "density", "force", "grav", "tend", "count"]
 
 SUBCOLOURS = {"none": "black",
               "density": "red",
               "force": "blue",
               "grav": "indigo",
+              "tend": "grey"
               "count": "purple"}
 
 #  Show docs if help is requested.

src/task.c

@@ -43,24 +43,28 @@
 /* Local headers. */
 #include "atomic.h"
 #include "error.h"
+#include "inline.h"
 #include "lock.h"
 
 /* Task type names. */
 const char *taskID_names[task_type_count] = {
-    "none",    "sort",          "self",          "pair",       "sub",
-    "init",    "ghost",         "drift",         "kick",       "kick_fixdt",
-    "send",    "recv",          "grav_gather_m", "grav_fft",   "grav_mm",
-    "grav_up", "grav_external", "part_sort",     "gpart_sort", "split_cell",
-    "rewait"};
+    "none",       "sort",    "self",          "pair",          "sub_self",
+    "sub_pair",   "init",    "ghost",         "drift",         "kick",
+    "kick_fixdt", "send",    "recv",          "grav_gather_m", "grav_fft",
+    "grav_mm",    "grav_up", "grav_external", "part_sort",     "gpart_sort",
+    "split_cell", "rewait"};
 
-const char *subtaskID_names[task_type_count] = {"none", "density", "force",
-                                                "grav"};
+const char *subtaskID_names[task_subtype_count] = {"none", "density", "force",
+                                                   "grav", "tend"};
 
 /**
  * @brief Computes the overlap between the parts array of two given cells.
  */
-size_t task_cell_overlap_part(const struct cell *ci, const struct cell *cj) {
+__attribute__((always_inline)) INLINE static size_t task_cell_overlap_part(
+    const struct cell *ci, const struct cell *cj) {
+
   if (ci == NULL || cj == NULL) return 0;
+
   if (ci->parts <= cj->parts &&
       ci->parts + ci->count >= cj->parts + cj->count) {
     return cj->count;
@@ -68,14 +72,18 @@ size_t task_cell_overlap_part(const struct cell *ci, const struct cell *cj) {
              cj->parts + cj->count >= ci->parts + ci->count) {
     return ci->count;
   }
+
   return 0;
 }
 
 /**
  * @brief Computes the overlap between the gparts array of two given cells.
  */
-size_t task_cell_overlap_gpart(const struct cell *ci, const struct cell *cj) {
+__attribute__((always_inline)) INLINE static size_t task_cell_overlap_gpart(
+    const struct cell *ci, const struct cell *cj) {
+
   if (ci == NULL || cj == NULL) return 0;
+
   if (ci->gparts <= cj->gparts &&
       ci->gparts + ci->gcount >= cj->gparts + cj->gcount) {
     return cj->gcount;
@@ -83,6 +91,7 @@ size_t task_cell_overlap_gpart(const struct cell *ci, const struct cell *cj) {
              cj->gparts + cj->gcount >= ci->gparts + ci->gcount) {
     return ci->gcount;
   }
+
   return 0;
 }
 
@@ -91,7 +100,8 @@ size_t task_cell_overlap_gpart(const struct cell *ci, const struct cell *cj) {
  *
  * @param t The #task.
  */
-enum task_actions task_acts_on(const struct task *t) {
+__attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
+    const struct task *t) {
 
   switch (t->type) {
 

src/task.h

@@ -112,6 +112,5 @@ int task_lock(struct task *t);
 void task_print_mask(unsigned int mask);
 void task_print_submask(unsigned int submask);
 void task_do_rewait(struct task *t);
-enum task_actions task_acts_on(const struct task *t);
 
 #endif /* SWIFT_TASK_H */

src/timestep.h

@@ -106,8 +106,9 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
 
     const float new_dt_external = gravity_compute_timestep_external(
         e->external_potential, e->physical_constants, p->gpart);
-    const float new_dt_self =
-        gravity_compute_timestep_self(e->physical_constants, p->gpart);
+    /* const float new_dt_self = */
+    /*     gravity_compute_timestep_self(e->physical_constants, p->gpart); */
+    const float new_dt_self = FLT_MAX; // MATTHIEU
 
     new_dt_grav = fminf(new_dt_external, new_dt_self);
   }

tests/test125cells.c

@@ -178,6 +178,21 @@ void get_solution(const struct cell *main_cell, struct solution_part *solution,
     solution[i].a_hydro[2] = -gradP[2] / solution[i].rho;
 
     solution[i].v_sig = 2.f * solution[i].c;
+
+    solution[i].S_dt = 0.f;
+    solution[i].u_dt = -(solution[i].P / solution[i].rho) * solution[i].div_v;
+  }
+}
+
+void reset_particles(struct cell *c, enum velocity_field vel,
+                     enum pressure_field press, float size, float density) {
+
+  for (size_t i = 0; i < c->count; ++i) {
+
+    struct part *p = &c->parts[i];
+
+    set_velocity(p, vel, size);
+    set_energy_state(p, press, size, density);
   }
 }
 
@@ -232,9 +247,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
         part->ti_begin = 0;
         part->ti_end = 1;
 
-        xpart->v_full[0] = part->v[0];
-        xpart->v_full[1] = part->v[1];
-        xpart->v_full[2] = part->v[2];
+        hydro_first_init_part(part, xpart);
         ++part;
         ++xpart;
       }
@@ -340,7 +353,8 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
               solution[pid].div_v, solution[pid].S, solution[pid].u,
               solution[pid].P, solution[pid].c, solution[pid].a_hydro[0],
               solution[pid].a_hydro[1], solution[pid].a_hydro[2],
-              solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt, 0.f);
+              solution[pid].h_dt, solution[pid].v_sig, solution[pid].S_dt,
+              solution[pid].u_dt);
     }
   }
 
@@ -427,6 +441,7 @@ int main(int argc, char *argv[]) {
 
   /* Help users... */
   message("Adiabatic index: ga = %f", hydro_gamma);
+  message("Hydro implementation: %s", SPH_IMPLEMENTATION);
   message("Smoothing length: h = %f", h * size);
   message("Kernel:               %s", kernel_name);
   message("Neighbour target: N = %f", h * h * h * kernel_norm);
@@ -447,17 +462,18 @@ int main(int argc, char *argv[]) {
   space.periodic = 0;
   space.h_max = h;
 
-  struct phys_const pc;
-  pc.const_newton_G = 1.f;
+  struct phys_const prog_const;
+  prog_const.const_newton_G = 1.f;
 
   struct hydro_props hp;
   hp.target_neighbours = h * h * h * kernel_norm;
   hp.delta_neighbours = 1.;
   hp.max_smoothing_iterations = 1;
+  hp.CFL_condition = 0.1;
 
   struct engine engine;
-  engine.physical_constants = &pc;
   engine.hydro_properties = &hp;
+  engine.physical_constants = &prog_const;
   engine.s = &space;
   engine.time = 0.1f;
   engine.ti_current = 1;
@@ -512,7 +528,7 @@ int main(int argc, char *argv[]) {
     for (int j = 0; j < 125; ++j) runner_do_init(&runner, cells[j], 0);
 
 /* Do the density calculation */
-#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
 
     /* Run all the pairs (only once !)*/
     for (int i = 0; i < 5; i++) {
@@ -554,7 +570,7 @@ int main(int argc, char *argv[]) {
     for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j]);
 
 /* Do the force calculation */
-#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
 
     /* Do the pairs (for the central 27 cells) */
     for (int i = 1; i < 4; i++) {
@@ -589,6 +605,9 @@ int main(int argc, char *argv[]) {
   /* Output timing */
   message("SWIFT calculation took       : %15lli ticks.", time / runs);
 
+  for (int j = 0; j < 125; ++j)
+    reset_particles(cells[j], vel, press, size, rho);
+
   /* NOW BRUTE-FORCE CALCULATION */
 
   const ticks tic = getticks();
@@ -597,7 +616,7 @@ int main(int argc, char *argv[]) {
   for (int j = 0; j < 125; ++j) runner_do_init(&runner, cells[j], 0);
 
 /* Do the density calculation */
-#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
 
   /* Run all the pairs (only once !)*/
   for (int i = 0; i < 5; i++) {
@@ -638,7 +657,7 @@ int main(int argc, char *argv[]) {
   for (int j = 0; j < 27; ++j) runner_do_ghost(&runner, inner_cells[j]);
 
 /* Do the force calculation */
-#if defined(DEFAULT_SPH) || !defined(WITH_VECTORIZATION)
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
 
   /* Do the pairs (for the central 27 cells) */
   for (int i = 1; i < 4; i++) {
@@ -654,6 +673,7 @@ int main(int argc, char *argv[]) {
 
   /* And now the self-interaction for the main cell */
   self_all_force(&runner, main_cell);
+
 #endif
 
   /* Finally, give a gentle kick */

tests/test125cells.sh.in

@@ -4,9 +4,18 @@ do
     for p in {0..2}
     do
 	echo ""
-	rm brute_force_125_standard.dat swift_dopair_125_standard.dat
+
+	rm -f brute_force_125_standard.dat swift_dopair_125_standard.dat
+
 	./test125cells -n 6 -r 1 -v $v -p $p -f standard
-	python @srcdir@/difffloat.py brute_force_125_standard.dat swift_dopair_125_standard.dat @srcdir@/tolerance_125.dat 6
+
+	if [ -e brute_force_125_standard.dat ]
+	then
+	    python @srcdir@/difffloat.py brute_force_125_standard.dat swift_dopair_125_standard.dat @srcdir@/tolerance_125.dat 6
+	else
+	    exit 1
+        fi
+
     done
 done
 	

tests/test27cells.c

@@ -318,6 +318,7 @@ int main(int argc, char *argv[]) {
 
   /* Help users... */
   message("Adiabatic index: ga = %f", hydro_gamma);
+  message("Hydro implementation: %s", SPH_IMPLEMENTATION);
   message("Smoothing length: h = %f", h * size);
   message("Kernel:               %s", kernel_name);
   message("Neighbour target: N = %f", h * h * h * kernel_norm);
@@ -366,7 +367,7 @@ int main(int argc, char *argv[]) {
 
     const ticks tic = getticks();
 
-#if defined(DEFAULT_SPH) || defined(GADGET2_SPH)
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
 
     /* Run all the pairs */
     for (int j = 0; j < 27; ++j)
@@ -402,7 +403,7 @@ int main(int argc, char *argv[]) {
 
   const ticks tic = getticks();
 
-#if defined(DEFAULT_SPH) || defined(GADGET2_SPH)
+#if !(defined(MINIMAL_SPH) && defined(WITH_VECTORIZATION))
 
   /* Run all the brute-force pairs */
   for (int j = 0; j < 27; ++j)

tests/test27cells.sh.in

@@ -3,9 +3,18 @@
 for v in {0..3}
 do
     echo ""
-    rm brute_force_27_standard.dat swift_dopair_27_standard.dat
+
+    rm -f brute_force_27_standard.dat swift_dopair_27_standard.dat
+
     ./test27cells -n 6 -r 1 -d 0 -f standard -v $v
-    python @srcdir@/difffloat.py brute_force_27_standard.dat swift_dopair_27_standard.dat @srcdir@/tolerance_27.dat 6
+
+    if [ -e brute_force_27_standard.dat ]
+    then
+	python @srcdir@/difffloat.py brute_force_27_standard.dat swift_dopair_27_standard.dat @srcdir@/tolerance_27_normal.dat 6
+    else
+	exit 1
+    fi
+
 done
 
 exit $?

tests/test27cellsPerturbed.sh.in

@@ -3,9 +3,18 @@
 for v in {0..3}
 do
     echo ""
-    rm brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat
+
+    rm -f brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat
+
     ./test27cells -n 6 -r 1 -d 0.1 -f perturbed -v $v
-    python @srcdir@/difffloat.py brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat @srcdir@/tolerance_27.dat 6
+
+    if [ -e brute_force_27_perturbed.dat ]
+    then
+	python @srcdir@/difffloat.py brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat @srcdir@/tolerance_27_perturbed.dat 6
+    else
+	exit 1
+    fi
+
 done
 
 exit $?

tests/testFFT.c

@@ -21,6 +21,12 @@
 #include <stdlib.h>
 #include <string.h>
 
+#ifndef HAVE_FFTW
+
+int main() { return 0; }
+
+#else
+
 #include <fftw3.h>
 
 /* Includes. */
@@ -187,3 +193,5 @@ int main() {
   free(gparts);
   return 0;
 }
+
+#endif

tests/testKernelGrav.c

@@ -74,10 +74,11 @@ int main() {
       swift_w = 1 / (r * r * r);
     }
 
-    printf("%2d: r= %f h= %f u= %f Wg(r,h)= %f Ws(r,h)= %f\n", k, r, h, u,
-           gadget_w, swift_w);
-
     if (fabsf(gadget_w - swift_w) > 2e-7) {
+
+      printf("%2d: r= %f h= %f u= %f Wg(r,h)= %f Ws(r,h)= %f\n", k, r, h, u,
+             gadget_w, swift_w);
+
       printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w);
       return 1;
     }
@@ -99,10 +100,11 @@ int main() {
 
     float gadget_w = erfc(u / 2) + u * exp(-u * u / 4) / sqrt(M_PI);
 
-    printf("%2d: r= %f r_lr= %f u= %f Ws(r)= %f Wg(r)= %f\n", k, r, a_smooth, u,
-           swift_w, gadget_w);
-
     if (fabsf(gadget_w - swift_w) > 2e-7) {
+
+      printf("%2d: r= %f r_lr= %f u= %f Ws(r)= %f Wg(r)= %f\n", k, r, a_smooth,
+             u, swift_w, gadget_w);
+
       printf("Invalid value ! Gadget= %e, SWIFT= %e\n", gadget_w, swift_w);
       return 1;
     }

tests/testMaths.c

@@ -39,24 +39,29 @@ int main() {
     const float rel =
         0.5f * fabs(exp_correct - exp_approx) / fabs(exp_correct + exp_approx);
 
-    printf("%2d: x= %f exp(x)= %e approx_exp(x)=%e abs=%e rel=%e\n", i, x,
-           exp_correct, exp_approx, abs, rel);
+    int error = 0;
 
     if (abs > 3e-6 && fabsf(x) <= 0.2) {
       printf("Absolute difference too large !\n");
-      return 1;
+      error = 1;
     }
     if (abs > 1.2e-7 && fabsf(x) <= 0.1) {
       printf("Absolute difference too large !\n");
-      return 1;
+      error = 1;
     }
 
     if (rel > 1e-6 && fabsf(x) <= 0.2) {
       printf("Relative difference too large !\n");
-      return 1;
+      error = 1;
     }
     if (rel > 4e-8 && fabsf(x) <= 0.1) {
       printf("Relative difference too large !\n");
+      error = 1;
+    }
+
+    if (error) {
+      printf("%2d: x= %f exp(x)= %e approx_exp(x)=%e abs=%e rel=%e\n", i, x,
+             exp_correct, exp_approx, abs, rel);
       return 1;
     }
   }

tests/testPair.sh.in

 #!/bin/bash
-rm brute_force_standard.dat swift_dopair_standard.dat
+
+echo ""
+
+rm -f brute_force_standard.dat swift_dopair_standard.dat
 
 ./testPair -p 6 -r 1 -d 0 -f standard
 
-python @srcdir@/difffloat.py brute_force_standard.dat swift_dopair_standard.dat @srcdir@/tolerance_pair.dat
+python @srcdir@/difffloat.py brute_force_standard.dat swift_dopair_standard.dat @srcdir@/tolerance_pair_normal.dat
 
 exit $?

tests/testPairPerturbed.sh.in

 #!/bin/bash
-rm brute_force_perturbed.dat swift_dopair_perturbed.dat
+
+echo ""
+
+rm -f brute_force_perturbed.dat swift_dopair_perturbed.dat
 
 ./testPair -p 6 -r 1 -d 0.1 -f perturbed
 
-python @srcdir@/difffloat.py brute_force_perturbed.dat swift_dopair_perturbed.dat @srcdir@/tolerance_pair.dat
+python @srcdir@/difffloat.py brute_force_perturbed.dat swift_dopair_perturbed.dat @srcdir@/tolerance_pair_perturbed.dat
 
 exit $?

tests/testParser.sh.in

 #!/bin/bash
 
-rm parser_output.yml
+rm -f parser_output.yml
 ./testParser @srcdir@/testParserInput.yaml

tests/tolerance_27_normal.dat

+#   ID      pos_x      pos_y      pos_z        v_x        v_y        v_z           rho        rho_dh        wcount     wcount_dh         div_v       curl_vx       curl_vy       curl_vz
+    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      2e-6	    2e-5       2e-3		 2e-6	     2e-6	   2e-6		 2e-6
+    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      1e-6	    1e-5       1e-4		 2e-5	     2e-5	   2e-5	 	 2e-5

tests/tolerance_27_perturbed.dat

+#   ID      pos_x      pos_y      pos_z        v_x        v_y        v_z           rho        rho_dh        wcount     wcount_dh         div_v       curl_vx       curl_vy       curl_vz
+    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1.2e-6     2e-6	    2.1e-5     2e-3		 2.1e-6	     2e-6	   2e-6		 2e-6
+    0	    1e-6       1e-6	  1e-6 	       1e-6 	  1e-6	     1e-6	   1e-6	      1e-6	    1e-5       1e-4		 2e-5	     4e-4	   4e-4	 	 4e-4