diff --git a/src/const.h b/src/const.h
index e6941f2cad4c62ed54f147628d47ccadab86050e..856043bf25f77e743de0daa99ea7665fc1534fe9 100644
--- a/src/const.h
+++ b/src/const.h
@@ -61,6 +61,6 @@
/* Gravity properties */
#define EXTERNAL_POTENTIAL_POINTMASS
-#define SANITY_CHECK
+#define SANITY_CHECKS
#endif /* SWIFT_CONST_H */
diff --git a/src/engine.c b/src/engine.c
index dd9f6cc5910a6049f2296e58b4917f7caecfdbd0..562e1aeed642bdcac3502073e302b4e8cb2f0b5c 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1066,6 +1066,69 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
#endif
}
+void engine_make_gravity_tasks(struct engine *e) {
+
+ struct space *s = e->s;
+ struct scheduler *sched = &e->sched;
+ const int nodeID = e->nodeID;
+ const int *cdim = s->cdim;
+ struct cell *cells = s->cells;
+
+ message("aa");
+
+ /* Run through the highest level of cells and add pairs. */
+ for (int i = 0; i < cdim[0]; i++) {
+ for (int j = 0; j < cdim[1]; j++) {
+ for (int k = 0; k < cdim[2]; k++) {
+
+ /* Get the cell */
+ const int cid = cell_getid(cdim, i, j, k);
+ struct cell *ci = &cells[cid];
+
+ /* Skip cells without gravity particles */
+ if (ci->gcount == 0) continue;
+
+ /* If the cells is local build a self-interaction */
+ if (ci->nodeID == nodeID)
+ scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci,
+ NULL, 0);
+
+ /* Now loop over all the neighbours of this cell */
+ for (int ii = -1; ii < 2; ii++) {
+ int iii = i + ii;
+ if (!s->periodic && (iii < 0 || iii >= cdim[0])) continue;
+ iii = (iii + cdim[0]) % cdim[0];
+ for (int jj = -1; jj < 2; jj++) {
+ int jjj = j + jj;
+ if (!s->periodic && (jjj < 0 || jjj >= cdim[1])) continue;
+ jjj = (jjj + cdim[1]) % cdim[1];
+ for (int kk = -1; kk < 2; kk++) {
+ int kkk = k + kk;
+ if (!s->periodic && (kkk < 0 || kkk >= cdim[2])) continue;
+ kkk = (kkk + cdim[2]) % cdim[2];
+
+ /* Get the neighbouring cell */
+ const int cjd = cell_getid(cdim, iii, jjj, kkk);
+ struct cell *cj = &cells[cjd];
+
+ /* Is that neighbour local and does it have particles ? */
+ if (cid >= cjd || cj->gcount == 0 ||
+ (ci->nodeID != nodeID && cj->nodeID != nodeID))
+ continue;
+
+ /* Construct the pair task */
+ const int sid =
+ sortlistID[(kk + 1) + 3 * ((jj + 1) + 3 * (ii + 1))];
+ scheduler_addtask(sched, task_type_pair, task_subtype_grav, sid,
+ 0, ci, cj, 1);
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
/**
* @brief Constructs the top-level pair tasks for the first hydro loop over
*neighbours
@@ -1435,8 +1498,7 @@ void engine_maketasks(struct engine *e) {
engine_make_hydroloop_tasks(e);
/* Add the gravity mm tasks. */
- /* if (e->policy & engine_policy_self_gravity) */
- /* engine_make_gravityinteraction_tasks(e); */
+ if (e->policy & engine_policy_self_gravity) engine_make_gravity_tasks(e);
/* Split the tasks. */
scheduler_splittasks(sched);
diff --git a/src/runner.c b/src/runner.c
index 5921691fba4c20adb136fee5513ae381761a6c71..1a094297aaa1fd66c347ae040fe43eaab4ba4e18 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -123,7 +123,6 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
if (g->ti_end <= ti_current) {
external_gravity(potential, constants, g);
-
}
}
if (timer) TIMER_TOC(timer_dograv_external);
@@ -1130,14 +1129,19 @@ void *runner_main(void *data) {
runner_doself1_density(r, ci);
else if (t->subtype == task_subtype_force)
runner_doself2_force(r, ci);
+ else if (t->subtype == task_subtype_grav)
+ runner_doself_grav(r, ci);
else
error("Unknown task subtype.");
break;
case task_type_pair:
+ message("bb");
if (t->subtype == task_subtype_density)
runner_dopair1_density(r, ci, cj);
else if (t->subtype == task_subtype_force)
runner_dopair2_force(r, ci, cj);
+ else if (t->subtype == task_subtype_grav)
+ runner_dopair_grav(r, ci, cj);
else
error("Unknown task subtype.");
break;
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 7e2293ba907ef6ae427f2b5e2ba5da156de12acf..113019a322df6ebf4c2cdefa24d44785f744c5f1 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -57,10 +57,10 @@ void runner_dograv_up(struct runner *r, struct cell *c) {
/**
* @brief Checks whether the cells are direct neighbours ot not. Both cells have
* to be of the same size
- *
+ *
*/
-__attribute__((always_inline)) INLINE static
-int are_neighbours(const struct cell *restrict ci, const struct cell *restrict cj) {
+__attribute__((always_inline)) INLINE static int are_neighbours(
+ const struct cell *restrict ci, const struct cell *restrict cj) {
#ifdef SANITY_CHECKS
if (ci->h[0] != cj->h[0])
@@ -80,8 +80,6 @@ int are_neighbours(const struct cell *restrict ci, const struct cell *restrict c
return 1;
}
-
-
/**
* @brief Computes the interaction of all the particles in a cell with the
* multipole of another cell.
@@ -90,8 +88,9 @@ int are_neighbours(const struct cell *restrict ci, const struct cell *restrict c
* @param ci The #cell with particles to interct.
* @param cj The #cell with the multipole.
*/
-__attribute__((always_inline)) INLINE static void runner_dograv_pair_pm(
- const struct runner *r, const struct cell *restrict ci, const struct cell *restrict cj) {
+__attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
+ const struct runner *r, const struct cell *restrict ci,
+ const struct cell *restrict cj) {
const struct engine *e = r->e;
const int gcount = ci->gcount;
@@ -101,7 +100,7 @@ __attribute__((always_inline)) INLINE static void runner_dograv_pair_pm(
TIMER_TIC;
-#ifdef SANITY_CHECK
+#ifdef SANITY_CHECKS
if (gcount == 0) error("Empty cell!"); // MATTHIEU sanity check
if (multi.mass == 0.0) // MATTHIEU sanity check
@@ -180,7 +179,7 @@ __attribute__((always_inline)) INLINE static void runner_dograv_pair_pm(
*
* @todo Use a local cache for the particles.
*/
-__attribute__((always_inline)) INLINE static void runner_dograv_pair_pp(
+__attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
struct runner *r, struct cell *ci, struct cell *cj) {
const struct engine *e = r->e;
@@ -192,7 +191,7 @@ __attribute__((always_inline)) INLINE static void runner_dograv_pair_pp(
TIMER_TIC;
-#ifdef SANITY_CHECK
+#ifdef SANITY_CHECKS
if (ci->h[0] != cj->h[0]) // MATTHIEU sanity check
error("Non matching cell sizes !! h_i=%f h_j=%f\n", ci->h[0], cj->h[0]);
#endif
@@ -250,7 +249,7 @@ __attribute__((always_inline)) INLINE static void runner_dograv_pair_pp(
* @todo Use a local cache for the particles.
*/
__attribute__((always_inline))
- INLINE static void runner_dograv_self_pp(struct runner *r, struct cell *c) {
+ INLINE static void runner_doself_grav_pp(struct runner *r, struct cell *c) {
const struct engine *e = r->e;
const int gcount = c->gcount;
@@ -259,7 +258,7 @@ __attribute__((always_inline))
TIMER_TIC;
-#ifdef SANITY_CHECK
+#ifdef SANITY_CHECKS
if (c->gcount == 0) // MATTHIEU sanity check
error("Empty cell !");
#endif
@@ -308,4 +307,109 @@ __attribute__((always_inline))
TIMER_TOC(TIMER_DOSELF); // MATTHIEU
}
+/**
+ * @brief Computes the interaction of all the particles in a cell with all the
+ * particles of another cell.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param cj The other #cell.
+ *
+ * @todo Use a local cache for the particles.
+ */
+static void runner_dopair_grav(struct runner *r, struct cell *ci,
+ struct cell *cj) {
+
+ const int gcount_i = ci->gcount;
+ const int gcount_j = cj->gcount;
+
+#ifdef SANITY_CHECKS
+
+ /* Early abort? */
+ if (gcount_i == 0 || gcount_j == 0) error("Empty cell !");
+
+ /* Bad stuff will happen if cell sizes are different */
+ if (ci->h[0] != cj->h[0])
+ error("Non matching cell sizes !! h_i=%f h_j=%f\n", ci->h[0], cj->h[0]);
+
+ /* Sanity check */
+ if (ci == cj)
+ error(
+ "The impossible has happened: pair interaction between a cell and "
+ "itself.");
+
+ /* Are the cells direct neighbours? */
+ if (!are_neighbours(ci, cj)) error("Non-neighbouring cells !");
+
+#endif
+
+ /* Are both cells split ? */
+ if (ci->split && cj->split) {
+
+ for (int j = 0; j < 8; j++) {
+ if (ci->progeny[j] != NULL) {
+
+ for (int k = 0; k < 8; k++) {
+ if (cj->progeny[k] != NULL) {
+
+ if (are_neighbours(ci->progeny[j], cj->progeny[k])) {
+
+ /* Recurse */
+ runner_dopair_grav(r, ci->progeny[j], cj->progeny[k]);
+
+ } else {
+
+ /* Ok, here we can go for particle-multipole interactions */
+ runner_dopair_grav_pm(r, ci->progeny[j], cj->progeny[k]);
+ runner_dopair_grav_pm(r, cj->progeny[k], ci->progeny[j]);
+ }
+ }
+ }
+ }
+ }
+ } else {/* Not split */
+
+ /* Compute the interactions at this level directly. */
+ runner_dopair_grav_pp(r, ci, cj);
+ }
+}
+
+static void runner_doself_grav(struct runner *r, struct cell *c) {
+
+ const int gcount = c->gcount;
+
+#ifdef SANITY_CHECKS
+
+ /* Early abort? */
+ if (gcount == 0) error("Empty cell !");
+#endif
+
+ message("aa");
+
+ /* If the cell is split, interact each progeny with itself, and with
+ each of its siblings. */
+ if (c->split) {
+
+ for (int j = 0; j < 8; j++) {
+ if (c->progeny[j] != NULL) {
+
+ runner_doself_grav(r, c->progeny[j]);
+
+ for (int k = j + 1; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+
+ runner_dopair_grav(r, c->progeny[j], c->progeny[k]);
+ }
+ }
+ }
+ }
+ }
+
+ /* If the cell is not split, then just go for it... */
+ else {
+
+ runner_doself_grav_pp(r, c);
+ }
+}
+
#endif /* SWIFT_RUNNER_DOIACT_GRAV_H */
diff --git a/src/scheduler.c b/src/scheduler.c
index 887e8d54fd1fda941fe0024c45733226ce8d4b43..539b36892adfc2a65d9d0c9954f06a2817eda9df 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -172,7 +172,8 @@ void scheduler_splittasks(struct scheduler *s) {
if (ci->split) {
/* Make a sub? */
- if (scheduler_dosub && ci->count < space_subsize / ci->count) {
+ if (scheduler_dosub && (ci->count * ci->count < space_subsize ||
+ ci->gcount * ci->gcount < space_subsize)) {
/* convert to a self-subtask. */
t->type = task_type_sub;
@@ -233,7 +234,8 @@ void scheduler_splittasks(struct scheduler *s) {
/* Replace by a single sub-task? */
if (scheduler_dosub &&
- ci->count * sid_scale[sid] < space_subsize / cj->count &&
+ (ci->count * cj->count * sid_scale[sid] < space_subsize ||
+ ci->gcount * cj->gcount * sid_scale[sid] < space_subsize) &&
sid != 0 && sid != 2 && sid != 6 && sid != 8) {
/* Make this task a sub task. */