diff --git a/src/engine.c b/src/engine.c
index b577ae3bf92234749802571f7b1d1c2462cb7d07..a21b5fe42d584f6e3ef77dfebdfbee1009d5e8ce 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1067,6 +1067,27 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
#endif
}
+__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])
+ error(" Cells of different size in distance calculation.");
+#endif
+
+ /* Maximum allowed distance */
+ const double min_dist = 1.2 * ci->h[0]; /* 1.2 accounts for rounding errors */
+
+ /* (Manhattan) Distance between the cells */
+ for (int k = 0; k < 3; k++) {
+ const double center_i = ci->loc[k];
+ const double center_j = cj->loc[k];
+ if (fabsf(center_i - center_j) > min_dist) return 0;
+ }
+
+ return 1;
+}
+
void engine_make_gravity_tasks(struct engine *e) {
struct space *s = e->s;
@@ -1101,8 +1122,20 @@ void engine_make_gravity_tasks(struct engine *e) {
/* Is that neighbour local ? */
if (cj->nodeID != nodeID) continue;
- scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci, cj,
- 1);
+#ifdef SANITY_CHECKS
+ if (ci == cj) {
+ message("oO");
+ continue;
+ }
+#endif
+
+ if (are_neighbours(ci, cj))
+ scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci,
+ cj, 1);
+ else
+ scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, 0, 0, ci,
+ cj, 1);
+
++counter;
}
}
@@ -1320,6 +1353,14 @@ void engine_make_gravity_dependencies(struct engine *e) {
/* Skip? */
if (t->skip) continue;
+ /* Long-range interaction */
+ if (t->type == task_type_grav_mm) {
+
+ scheduler_addunlock(sched, t->ci->super->init, t);
+ scheduler_addunlock(sched, t->ci->super->grav_up, t);
+ scheduler_addunlock(sched, t, t->ci->super->kick);
+ }
+
/* Self-interaction? */
if (t->type == task_type_self && t->subtype == task_subtype_grav) {
@@ -2296,6 +2337,7 @@ void engine_step(struct engine *e) {
if (e->policy & engine_policy_self_gravity) {
mask |= 1 << task_type_grav_up;
+ mask |= 1 << task_type_grav_mm;
mask |= 1 << task_type_self;
mask |= 1 << task_type_pair;
mask |= 1 << task_type_sub;
@@ -2331,6 +2373,9 @@ void engine_step(struct engine *e) {
s->gparts[k].x[0], s->gparts[k].x[1], s->gparts[k].x[2],
s->gparts[k].a_grav[0], s->gparts[k].a_grav[1],
s->gparts[k].a_grav[2]);
+ if (s->gparts[k].id == -1)
+ message("interacting mass= %f (%f)", s->gparts[k].mass_interacted,
+ s->gparts[k].mass_interacted + s->gparts[k].mass);
}
fclose(file);
@@ -2339,10 +2384,10 @@ void engine_step(struct engine *e) {
memcpy(temp, s->gparts, s->nr_gparts * sizeof(struct gpart));
gravity_n2(temp, s->nr_gparts);
file = fopen("grav_brute.dat", "w");
- for(size_t k = 0; k < s->nr_gparts; ++k) {
- fprintf(file, "%lld %f %f %f %f %f %f\n", temp[k].id,
- temp[k].x[0], temp[k].x[1], temp[k].x[2],
- temp[k].a_grav[0], temp[k].a_grav[1], temp[k].a_grav[2]);
+ for (size_t k = 0; k < s->nr_gparts; ++k) {
+ fprintf(file, "%lld %f %f %f %f %f %f\n", temp[k].id, temp[k].x[0],
+ temp[k].x[1], temp[k].x[2], temp[k].a_grav[0], temp[k].a_grav[1],
+ temp[k].a_grav[2]);
}
fclose(file);
diff --git a/src/runner.c b/src/runner.c
index 4d32cfd24def855858eb493382aefaa09ee84ac5..49223b46a11c1b113b9f0633ccc961e70b9bee43 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1188,11 +1188,11 @@ void *runner_main(void *data) {
/* else */
/* runner_dopair_grav(r, t->ci, t->cj); */
/* break; */
- /* case task_type_grav_mm: */
- /* runner_dograv_mm(r, t->ci, t->cj); */
- /* break; */
+ case task_type_grav_mm:
+ runner_do_grav_mm(r, t->ci, t->cj);
+ break;
case task_type_grav_up:
- runner_dograv_up(r, t->ci);
+ runner_do_grav_up(r, t->ci);
break;
/* case task_type_grav_down: */
/* runner_dograv_down(r, t->ci); */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 9ec64f74f1f41e8a68a36e1efa42f7575867a4c5..2bed71825da5aa608ccfdc8f8ebd6ab2af7d3366 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -25,7 +25,7 @@
#include "clocks.h"
#include "part.h"
-#define ICHECK 1
+#define ICHECK -1
/**
* @brief Compute the recursive upward sweep, i.e. construct the
@@ -34,13 +34,13 @@
* @param r The #runner.
* @param c The top-level #cell.
*/
-void runner_dograv_up(struct runner *r, struct cell *c) {
+void runner_do_grav_up(struct runner *r, struct cell *c) {
if (c->split) {/* Regular node */
/* Recurse. */
for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL) runner_dograv_up(r, c->progeny[k]);
+ if (c->progeny[k] != NULL) runner_do_grav_up(r, c->progeny[k]);
/* Collect the multipoles from the progeny. */
multipole_reset(&c->multipole);
@@ -98,11 +98,11 @@ __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 struct engine *e = r->e;
const int gcount = ci->gcount;
struct gpart *restrict gparts = ci->gparts;
const struct multipole multi = cj->multipole;
- const int ti_current = e->ti_current;
+ // const int ti_current = e->ti_current;
TIMER_TIC;
@@ -115,6 +115,17 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
/* Anything to do here? */
// if (ci->ti_end_min > ti_current) return;
+ /* Loop over every particle in leaf. */
+ /* for (int pid = 0; pid < gcount; pid++) { */
+
+ /* /\* Get a hold of the ith part in ci. *\/ */
+ /* struct gpart *restrict gp = &gparts[pid]; */
+
+ /* if (gp->id == -ICHECK) */
+ /* message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id,
+ * cj->loc[0], */
+ /* cj->loc[1], cj->loc[2], cj->h[0], cj->gcount); */
+ /* } */
/* Loop over every particle in leaf. */
for (int pid = 0; pid < gcount; pid++) {
@@ -122,12 +133,9 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
/* Get a hold of the ith part in ci. */
struct gpart *restrict gp = &gparts[pid];
- if (gp->ti_end > ti_current) continue;
+ if (gp->id == -ICHECK) message("id=%lld mass= %f", gp->id, multi.mass);
- if (gp->id == -ICHECK)
- message("id=%lld x=[%f %f %f] mass=%f CoM=[ %f %f %f ] size= %f\n", gp->id,
- gp->x[0], gp->x[1], gp->x[2], multi.mass,
- multi.CoM[0], multi.CoM[1], multi.CoM[2], cj->h[0]);
+ // if (gp->ti_end > ti_current) continue;
/* Compute the pairwise distance. */
const float dx[3] = {multi.CoM[0] - gp->x[0], // x
@@ -147,7 +155,6 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
gp->a_grav[2] += mrinv3 * dx[2];
#elif multipole_order == 2
-
/* Terms up to 2nd order (quadrupole) */
/* Follows the notation in Bonsai */
@@ -194,18 +201,18 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
__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;
+ // const struct engine *e = r->e;
const int gcount_i = ci->gcount;
const int gcount_j = cj->gcount;
struct gpart *restrict gparts_i = ci->gparts;
struct gpart *restrict gparts_j = cj->gparts;
- const int ti_current = e->ti_current;
+ // const int ti_current = e->ti_current;
TIMER_TIC;
#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]);
+ error("Non matching cell sizes !! h_i=%f h_j=%f", ci->h[0], cj->h[0]);
#endif
/* Anything to do here? */
@@ -213,21 +220,24 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
for (int pid = 0; pid < gcount_i; pid++) {
+ /* Get a hold of the ith part in ci. */
struct gpart *restrict gp = &gparts_i[pid];
- if (gp->id == -ICHECK)
- message("id=%lld size=%f\n", gp->id, cj->h[0]);
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, cj->loc[0],
+ cj->loc[1], cj->loc[2], cj->h[0], cj->gcount);
}
for (int pid = 0; pid < gcount_j; pid++) {
+ /* Get a hold of the ith part in ci. */
struct gpart *restrict gp = &gparts_j[pid];
- if (gp->id == -ICHECK)
- message("id=%lld size=%f\n", gp->id, ci->h[0]);
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count=%d", gp->id, ci->loc[0],
+ ci->loc[1], ci->loc[2], ci->h[0], ci->gcount);
}
-
/* Loop over all particles in ci... */
for (int pid = 0; pid < gcount_i; pid++) {
@@ -253,18 +263,18 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
const float w = ir * ir * ir;
const float wdx[3] = {w * dx[0], w * dx[1], w * dx[2]};
- if (gpi->ti_end <= ti_current) {
- gpi->a_grav[0] -= wdx[0] * mj;
- gpi->a_grav[1] -= wdx[1] * mj;
- gpi->a_grav[2] -= wdx[2] * mj;
- gpi->mass_interacted += mj;
- }
- if (gpj->ti_end <= ti_current) {
- gpj->a_grav[0] += wdx[0] * mi;
- gpj->a_grav[1] += wdx[1] * mi;
- gpj->a_grav[2] += wdx[2] * mi;
- gpj->mass_interacted += mi;
- }
+ // if (gpi->ti_end <= ti_current) {
+ gpi->a_grav[0] -= wdx[0] * mj;
+ gpi->a_grav[1] -= wdx[1] * mj;
+ gpi->a_grav[2] -= wdx[2] * mj;
+ gpi->mass_interacted += mj;
+ //}
+ // if (gpj->ti_end <= ti_current) {
+ gpj->a_grav[0] += wdx[0] * mi;
+ gpj->a_grav[1] += wdx[1] * mi;
+ gpj->a_grav[2] += wdx[2] * mi;
+ gpj->mass_interacted += mi;
+ // }
}
}
@@ -282,10 +292,10 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
__attribute__((always_inline))
INLINE static void runner_doself_grav_pp(struct runner *r, struct cell *c) {
- const struct engine *e = r->e;
+ // const struct engine *e = r->e;
const int gcount = c->gcount;
struct gpart *restrict gparts = c->gparts;
- const int ti_current = e->ti_current;
+ // const int ti_current = e->ti_current;
TIMER_TIC;
@@ -297,6 +307,16 @@ __attribute__((always_inline))
/* Anything to do here? */
// if (c->ti_end_min > ti_current) return;
+ for (int pid = 0; pid < gcount; pid++) {
+
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &gparts[pid];
+
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, c->loc[0],
+ c->loc[1], c->loc[2], c->h[0], c->gcount);
+ }
+
/* Loop over all particles in ci... */
for (int pid = 0; pid < gcount; pid++) {
@@ -322,18 +342,18 @@ __attribute__((always_inline))
const float w = ir * ir * ir;
const float wdx[3] = {w * dx[0], w * dx[1], w * dx[2]};
- if (gpi->ti_end <= ti_current) {
- gpi->a_grav[0] -= wdx[0] * mj;
- gpi->a_grav[1] -= wdx[1] * mj;
- gpi->a_grav[2] -= wdx[2] * mj;
- gpi->mass_interacted += mj;
- }
- if (gpj->ti_end <= ti_current) {
- gpj->a_grav[0] += wdx[0] * mi;
- gpj->a_grav[1] += wdx[1] * mi;
- gpj->a_grav[2] += wdx[2] * mi;
- gpj->mass_interacted += mi;
- }
+ // if (gpi->ti_end <= ti_current) {
+ gpi->a_grav[0] -= wdx[0] * mj;
+ gpi->a_grav[1] -= wdx[1] * mj;
+ gpi->a_grav[2] -= wdx[2] * mj;
+ gpi->mass_interacted += mj;
+ //}
+ // if (gpj->ti_end <= ti_current) {
+ gpj->a_grav[0] += wdx[0] * mi;
+ gpj->a_grav[1] += wdx[1] * mi;
+ gpj->a_grav[2] += wdx[2] * mi;
+ gpj->mass_interacted += mi;
+ //}
}
}
@@ -363,7 +383,7 @@ static void runner_dopair_grav(struct runner *r, struct cell *ci,
/* 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]);
+ error("Non matching cell sizes !! h_i=%f h_j=%f", ci->h[0], cj->h[0]);
/* Sanity check */
if (ci == cj)
@@ -373,68 +393,69 @@ static void runner_dopair_grav(struct runner *r, struct cell *ci,
#endif
- /* Are the cells direct neighbours? */
- if (!are_neighbours(ci, cj)) {
+ for (int pid = 0; pid < gcount_i; pid++) {
- /* Ok, here we can go for particle-multipole interactions */
- runner_dopair_grav_pm(r, ci, cj);
- runner_dopair_grav_pm(r, cj, ci);
-
- /* error( */
- /* "Non-neighbouring cells ! ci->x=[%f %f %f] ci->h=%f cj->loc=[%f %f
- * %f] " */
- /* "cj->h=%f", */
- /* ci->loc[0], ci->loc[1], ci->loc[2], ci->h[0], cj->loc[0], cj->loc[1],
- */
- /* cj->loc[2], cj->h[0]); */
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &ci->gparts[pid];
+
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, cj->loc[0],
+ cj->loc[1], cj->loc[2], cj->h[0], cj->gcount);
}
- /* for (int i = 0; i < gcount_i; i++) { */
+ for (int pid = 0; pid < gcount_j; pid++) {
- /* struct gpart *const gp = &ci->gparts[i]; */
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &cj->gparts[pid];
- /* if (gp->id == -ICHECK) */
- /* message("id=%lld a=[%f %f %f]\n", gp->id, gp->a_grav[0], gp->a_grav[1], */
- /* gp->a_grav[2]); */
- /* } */
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, ci->loc[0],
+ ci->loc[1], ci->loc[2], ci->h[0], ci->gcount);
+ }
- /* for (int i = 0; i < gcount_j; i++) { */
+ /* Are the cells direct neighbours? */
+ if (!are_neighbours(ci, cj)) {
- /* struct gpart *const gp = &cj->gparts[i]; */
+ /* Ok, here we can go for particle-multipole interactions */
+ // runner_dopair_grav_pm(r, ci, cj);
+ // runner_dopair_grav_pm(r, cj, ci);
- /* if (gp->id == -ICHECK) */
- /* message("id=%lld a=[%f %f %f]\n", gp->id, gp->a_grav[0], gp->a_grav[1], */
- /* gp->a_grav[2]); */
- /* } */
+ error(
+ "Non-neighbouring cells ! ci->x=[%f %f %f] ci->h=%f cj->loc=[%f %f %f] "
+ "cj->h=%f",
+ ci->loc[0], ci->loc[1], ci->loc[2], ci->h[0], cj->loc[0], cj->loc[1],
+ cj->loc[2], cj->h[0]);
+ } else {
- /* Are both cells split ? */
- if (ci->split && cj->split) {
+ /* Are both cells split ? */
+ if (ci->split && cj->split) {
- for (int j = 0; j < 8; j++) {
- if (ci->progeny[j] != NULL) {
+ for (int j = 0; j < 8; j++) {
+ if (ci->progeny[j] != NULL) {
- for (int k = 0; k < 8; k++) {
- if (cj->progeny[k] != NULL) {
+ for (int k = 0; k < 8; k++) {
+ if (cj->progeny[k] != NULL) {
- if (are_neighbours(ci->progeny[j], cj->progeny[k])) {
+ if (are_neighbours(ci->progeny[j], cj->progeny[k])) {
- /* Recurse */
- runner_dopair_grav(r, ci->progeny[j], cj->progeny[k]);
+ /* Recurse */
+ runner_dopair_grav(r, ci->progeny[j], cj->progeny[k]);
- } else {
+ } 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]);
+ /* 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 */
+ } else {/* Not split */
- /* Compute the interactions at this level directly. */
- runner_dopair_grav_pp(r, ci, cj);
+ /* Compute the interactions at this level directly. */
+ runner_dopair_grav_pp(r, ci, cj);
+ }
}
}
@@ -448,14 +469,15 @@ static void runner_doself_grav(struct runner *r, struct cell *c) {
if (gcount == 0) error("Empty cell !");
#endif
- for (int i = 0; i < gcount; i++) {
+ /* for (int i = 0; i < gcount; i++) { */
- struct gpart *const gp = &c->gparts[i];
+ /* struct gpart *const gp = &c->gparts[i]; */
- if (gp->id == -ICHECK)
- message("id=%lld a=[%f %f %f]\n", gp->id, gp->a_grav[0], gp->a_grav[1],
- gp->a_grav[2]);
- }
+ /* if (gp->id == -ICHECK) */
+ /* message("id=%lld a=[%f %f %f]", gp->id, gp->a_grav[0], gp->a_grav[1],
+ */
+ /* gp->a_grav[2]); */
+ /* } */
/* If the cell is split, interact each progeny with itself, and with
each of its siblings. */
@@ -493,8 +515,10 @@ static void runner_dosub_grav(struct runner *r, struct cell *ci,
} else {
- /* if (!are_neighbours(ci, cj)) { */
+ if (!are_neighbours(ci, cj)) {
+ error("Non-neighbouring cells in pair task !");
+ }
/* /\* Ok, here we can go for particle-multipole interactions *\/ */
/* runner_dopair_grav_pm(r, ci, cj); */
/* runner_dopair_grav_pm(r, cj, ci); */
@@ -508,4 +532,36 @@ static void runner_dosub_grav(struct runner *r, struct cell *ci,
}
}
+static void runner_do_grav_mm(struct runner *r, struct cell *ci,
+ struct cell *cj) {
+
+ for (int pid = 0; pid < ci->gcount; pid++) {
+
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &ci->gparts[pid];
+
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, cj->loc[0],
+ cj->loc[1], cj->loc[2], cj->h[0], cj->gcount);
+ }
+
+ for (int pid = 0; pid < cj->gcount; pid++) {
+
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gp = &cj->gparts[pid];
+
+ if (gp->id == -ICHECK)
+ message("id=%lld loc=[ %f %f %f ] size= %f count= %d", gp->id, ci->loc[0],
+ ci->loc[1], ci->loc[2], ci->h[0], ci->gcount);
+ }
+
+ if (are_neighbours(ci, cj)) {
+
+ error("Non-neighbouring cells in mm task !");
+ }
+
+ runner_dopair_grav_pm(r, ci, cj);
+ runner_dopair_grav_pm(r, cj, ci);
+}
+
#endif /* SWIFT_RUNNER_DOIACT_GRAV_H */
diff --git a/src/scheduler.c b/src/scheduler.c
index 6271f80ca306f3546f12360f18eb7db038d97fa3..44859d4c52de1b2b43f1258526b14b646bdf1fe8 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -172,14 +172,12 @@ void scheduler_splittasks(struct scheduler *s) {
if (ci->split) {
/* Make a sub? */
- if (scheduler_dosub && (ci->count * ci->count < space_subsize)) {
- // ci->gcount * ci->gcount < space_subsize)) {
+ 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;
- // message("TASK type=%s subtype=%s", taskID_names[t->type],
- // subtaskID_names[t->subtype]);
}
/* Otherwise, make tasks explicitly. */
@@ -188,8 +186,6 @@ void scheduler_splittasks(struct scheduler *s) {
/* Take a step back (we're going to recycle the current task)... */
redo = 1;
- // message("aa");
-
/* Add the self task. */
int first_child = 0;
while (ci->progeny[first_child] == NULL) first_child++;
@@ -210,9 +206,39 @@ void scheduler_splittasks(struct scheduler *s) {
}
}
+ /* /\* Hydro Pair interaction? *\/ */
+ /* else if (t->type == task_type_pair && t->subtype == task_subtype_grav) {
+ */
+
+ /* /\* Get a handle on the cells involved. *\/ */
+ /* struct cell *ci = t->ci; */
+ /* struct cell *cj = t->cj; */
+
+ /* /\* Foreign task? *\/ */
+ /* if (ci->nodeID != s->nodeID && cj->nodeID != s->nodeID) { */
+ /* t->skip = 1; */
+ /* continue; */
+ /* } */
+
+ /* if (ci->split && cj->split ) { */
+
+ /* /\* Replace by a single sub-task? *\/ */
+ /* if (scheduler_dosub && ci->gcount * cj->gcount < space_subsize) { */
+
+ /* /\* Make this task a sub task. *\/ */
+ /* t->type = task_type_sub; */
+ /* } */
+
+ /* else { */
+ /* message("oO"); */
+ /* } */
- /* Pair interaction? */
- else if (t->type == task_type_pair) {
+ /* } */
+
+ /* } */
+
+ /* Hydro Pair interaction? */
+ else if (t->type == task_type_pair && t->subtype != task_subtype_grav) {
/* Get a handle on the cells involved. */
struct cell *ci = t->ci;
@@ -238,8 +264,7 @@ void scheduler_splittasks(struct scheduler *s) {
/* Replace by a single sub-task? */
if (scheduler_dosub &&
- (ci->count * cj->count * sid_scale[sid] < space_subsize) &&// ||
- //ci->gcount * cj->gcount * sid_scale[sid] < space_subsize) &&
+ ci->count * cj->count * sid_scale[sid] < space_subsize &&
sid != 0 && sid != 2 && sid != 6 && sid != 8) {
/* Make this task a sub task. */
@@ -523,139 +548,17 @@ void scheduler_splittasks(struct scheduler *s) {
} /* pair interaction? */
- /* Gravity interaction? */
- /* else if (t->type == task_type_grav_mm) { */
+ /* Long-range gravity interaction ? */
+ else if (t->type == task_type_grav_mm) {
- /* /\* Get a handle on the cells involved. *\/ */
- /* struct cell *ci = t->ci; */
- /* struct cell *cj = t->cj; */
-
- /* /\* Self-interaction? *\/ */
- /* if (cj == NULL) { */
-
- /* /\* Ignore this task if the cell has no gparts. *\/ */
- /* if (ci->gcount == 0) t->type = task_type_none; */
-
- /* /\* If the cell is split, recurse. *\/ */
- /* else if (ci->split) { */
-
- /* /\* Make a single sub-task? *\/ */
- /* if (scheduler_dosub && ci->gcount < space_subsize / ci->gcount) {
- */
-
- /* t->type = task_type_sub; */
- /* t->subtype = task_subtype_grav; */
-
- /* } */
-
- /* /\* Otherwise, just split the task. *\/ */
- /* else { */
-
- /* /\* Split this task into tasks on its progeny. *\/ */
- /* t->type = task_type_none; */
- /* for (int j = 0; j < 8; j++) */
- /* if (ci->progeny[j] != NULL && ci->progeny[j]->gcount > 0) { */
- /* if (t->type == task_type_none) { */
- /* t->type = task_type_grav_mm; */
- /* t->ci = ci->progeny[j]; */
- /* t->cj = NULL; */
- /* } else */
- /* t = scheduler_addtask(s, task_type_grav_mm,
- * task_subtype_none, */
- /* 0, 0, ci->progeny[j], NULL, 0); */
- /* for (int k = j + 1; k < 8; k++) */
- /* if (ci->progeny[k] != NULL && ci->progeny[k]->gcount > 0) {
- */
- /* if (t->type == task_type_none) { */
- /* t->type = task_type_grav_mm; */
- /* t->ci = ci->progeny[j]; */
- /* t->cj = ci->progeny[k]; */
- /* } else */
- /* t = scheduler_addtask(s, task_type_grav_mm, */
- /* task_subtype_none, 0, 0, */
- /* ci->progeny[j], ci->progeny[k],
- * 0); */
- /* } */
- /* } */
- /* redo = (t->type != task_type_none); */
- /* } */
-
- /* } */
-
- /* /\* Otherwise, just make a pp task out of it. *\/ */
- /* else */
- /* t->type = task_type_grav_pp; */
+ /* Get a handle on the cells involved. */
+ struct cell *ci = t->ci;
+ struct cell *cj = t->cj;
- /* } */
+ /* Safety thing */
+ if (ci->gcount == 0 || cj->gcount == 0) t->type = task_type_none;
- /* /\* Nope, pair. *\/ */
- /* else { */
-
- /* /\* Make a sub-task? *\/ */
- /* if (scheduler_dosub && ci->gcount < space_subsize / cj->gcount) { */
-
- /* t->type = task_type_sub; */
- /* t->subtype = task_subtype_grav; */
-
- /* } */
-
- /* /\* Otherwise, split the task. *\/ */
- /* else { */
-
- /* /\* Get the opening angle theta. *\/ */
- /* float dx[3], theta; */
- /* for (int k = 0; k < 3; k++) { */
- /* dx[k] = fabs(ci->loc[k] - cj->loc[k]); */
- /* if (s->space->periodic && dx[k] > 0.5 * s->space->dim[k]) */
- /* dx[k] = -dx[k] + s->space->dim[k]; */
- /* if (dx[k] > 0.0f) dx[k] -= ci->h[k]; */
- /* } */
- /* theta = */
- /* (dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]) / */
- /* (ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] *
- * ci->h[2]); */
-
- /* /\* Ignore this task if the cell has no gparts. *\/ */
- /* if (ci->gcount == 0 || cj->gcount == 0) t->type = task_type_none;
- */
-
- /* /\* Split the interaction? *\/ */
- /* else if (theta < const_theta_max * const_theta_max) { */
-
- /* /\* Are both ci and cj split? *\/ */
- /* if (ci->split && cj->split) { */
-
- /* /\* Split this task into tasks on its progeny. *\/ */
- /* t->type = task_type_none; */
- /* for (int j = 0; j < 8; j++) */
- /* if (ci->progeny[j] != NULL && ci->progeny[j]->gcount > 0) {
- */
- /* for (int k = 0; k < 8; k++) */
- /* if (cj->progeny[k] != NULL && cj->progeny[k]->gcount > 0)
- * { */
- /* if (t->type == task_type_none) { */
- /* t->type = task_type_grav_mm; */
- /* t->ci = ci->progeny[j]; */
- /* t->cj = cj->progeny[k]; */
- /* } else */
- /* t = scheduler_addtask( */
- /* s, task_type_grav_mm, task_subtype_none, 0, 0, */
- /* ci->progeny[j], cj->progeny[k], 0); */
- /* } */
- /* } */
- /* redo = (t->type != task_type_none); */
-
- /* } */
-
- /* /\* Otherwise, make a pp task out of it. *\/ */
- /* else */
- /* t->type = task_type_grav_pp; */
- /* } */
- /* } */
-
- /* } /\* gravity pair interaction? *\/ */
-
- /* } /\* gravity interaction? *\/ */
+ } /* gravity interaction? */
} /* loop over all tasks. */
}
diff --git a/src/space.c b/src/space.c
index d87107a80add1b04674b74a8af55004427ddeac3..c4b81646f8d8e40504e5d18332ee547af0a7fefc 100644
--- a/src/space.c
+++ b/src/space.c
@@ -180,6 +180,8 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
}
}
+ message("h_max: %f", s->h_max);
+
/* If we are running in parallel, make sure everybody agrees on
how large the largest cell should be. */
#ifdef WITH_MPI
diff --git a/src/task.c b/src/task.c
index 77503d2207ae1a80e1c16f9eecc9fd59ed93f871..4069ac520f1864ccb230b517d05281e42d491c61 100644
--- a/src/task.c
+++ b/src/task.c
@@ -47,10 +47,10 @@
/* Task type names. */
const char *taskID_names[task_type_count] = {
- "none", "sort", "self", "pair", "sub",
- "init", "ghost", "drift", "kick", "send",
- "recv", "grav_up", "grav_external", "part_sort", "gpart_sort",
- "split_cell", "rewait"};
+ "none", "sort", "self", "pair", "sub",
+ "init", "ghost", "drift", "kick", "send",
+ "recv", "grav_mm", "grav_up", "grav_external", "part_sort",
+ "gpart_sort", "split_cell", "rewait"};
const char *subtaskID_names[task_type_count] = {"none", "density",
"force", "grav"};
diff --git a/src/task.h b/src/task.h
index 4956860c96020a211db5137c96ad4d260582e1c7..b159fe0e578695a0f20cbe0744397a7b0d4d3bef 100644
--- a/src/task.h
+++ b/src/task.h
@@ -45,7 +45,7 @@ enum task_types {
task_type_send,
task_type_recv,
/* task_type_grav_pp, */
- /* task_type_grav_mm, */
+ task_type_grav_mm,
task_type_grav_up,
/* task_type_grav_down, */
task_type_grav_external,