diff --git a/examples/UniformDMBox/makeIC.py b/examples/UniformDMBox/makeIC.py
index 449d780fb31bc23dd194f772be45d35e6b0bbe3f..2aee89798a5b8bbd425a6b73528779fb1aa7db23 100644
--- a/examples/UniformDMBox/makeIC.py
+++ b/examples/UniformDMBox/makeIC.py
@@ -26,7 +26,7 @@ from numpy import *
# with a density of 1
# Parameters
-periodic= 1 # 1 For periodic box
+periodic= 0 # 1 For periodic box
boxSize = 1.
rho = 1.
L = int(sys.argv[1]) # Number of particles along one axis
diff --git a/src/engine.c b/src/engine.c
index 562e1aeed642bdcac3502073e302b4e8cb2f0b5c..f7e749ab109a9223da7af7020b7b34ad0e5c6278 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -62,6 +62,7 @@
#include "serial_io.h"
#include "single_io.h"
#include "timers.h"
+#include "tools.h"
const char *engine_policy_names[13] = {
"none", "rand", "steal", "keep",
@@ -1271,6 +1272,68 @@ void engine_count_and_link_tasks(struct engine *e) {
}
}
+void engine_make_gravity_dependencies(struct engine *e) {
+
+ struct scheduler *sched = &e->sched;
+ const int nodeID = e->nodeID;
+ const int nr_tasks = sched->nr_tasks;
+
+ for (int k = 0; k < nr_tasks; k++) {
+
+ /* Get a pointer to the task. */
+ struct task *t = &sched->tasks[k];
+
+ /* Skip? */
+ if (t->skip) continue;
+
+ /* Self-interaction? */
+ if (t->type == task_type_self && t->subtype == task_subtype_grav) {
+
+ 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);
+
+ }
+
+ /* Otherwise, pair interaction? */
+ else if (t->type == task_type_pair && t->subtype == task_subtype_grav) {
+
+ if (t->ci->nodeID == nodeID) {
+
+ 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);
+ }
+
+ if (t->cj->nodeID == nodeID && t->ci->super != t->cj->super) {
+
+ scheduler_addunlock(sched, t->cj->super->init, t);
+ scheduler_addunlock(sched, t->cj->super->grav_up, t);
+ scheduler_addunlock(sched, t, t->cj->super->kick);
+ }
+
+ }
+
+ /* Otherwise, sub interaction? */
+ else if (t->type == task_type_sub && t->subtype == task_subtype_grav) {
+
+ if (t->ci->nodeID == nodeID) {
+
+ 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);
+ }
+ if (t->cj != NULL && t->cj->nodeID == nodeID &&
+ t->ci->super != t->cj->super) {
+
+ scheduler_addunlock(sched, t->cj->super->init, t);
+ scheduler_addunlock(sched, t->cj->super->grav_up, t);
+ scheduler_addunlock(sched, t, t->cj->super->kick);
+ }
+ }
+ }
+}
+
/**
* @brief Creates the dependency network for the hydro tasks of a given cell.
*
@@ -1382,10 +1445,6 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
}
}
- /* /\* Kick tasks should rely on the grav_down tasks of their cell. *\/ */
- /* else if (t->type == task_type_kick && t->ci->grav_down != NULL) */
- /* scheduler_addunlock(sched, t->ci->grav_down, t); */
-
/* External gravity tasks should depend on init and unlock the kick */
else if (t->type == task_type_grav_external) {
scheduler_addunlock(sched, t->ci->init, t);
@@ -1530,6 +1589,10 @@ void engine_maketasks(struct engine *e) {
of its super-cell. */
engine_make_extra_hydroloop_tasks(e);
+ /* Add the dependencies for the self-gravity stuff */
+ if (e->policy & engine_policy_self_gravity)
+ engine_make_gravity_dependencies(e);
+
#ifdef WITH_MPI
/* Add the communication tasks if MPI is being used. */
@@ -2011,6 +2074,10 @@ void engine_init_particles(struct engine *e) {
if (e->policy & engine_policy_self_gravity) {
mask |= 1 << task_type_grav_up;
+ mask |= 1 << task_type_self;
+ mask |= 1 << task_type_pair;
+
+ submask |= 1 << task_subtype_grav;
}
/* Add the tasks corresponding to external gravity to the masks */
@@ -2192,6 +2259,10 @@ void engine_step(struct engine *e) {
if (e->policy & engine_policy_self_gravity) {
mask |= 1 << task_type_grav_up;
+ mask |= 1 << task_type_self;
+ mask |= 1 << task_type_pair;
+
+ submask |= 1 << task_subtype_grav;
}
/* Add the tasks corresponding to external gravity to the masks */
@@ -2216,6 +2287,31 @@ void engine_step(struct engine *e) {
/* Check that the multipoles are correct */
space_map_cells_pre(s, 1, cell_check_multipole, NULL);
+ FILE *file = fopen("grav_swift.dat", "w");
+ for (size_t k = 0; k < s->nr_gparts; ++k) {
+ fprintf(file, "%lld %f %f %f %f %f %f\n", s->gparts[k].id,
+ 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]);
+ }
+ fclose(file);
+
+ /* Check the gravity accelerations */
+ /* struct gpart *temp = malloc(s->nr_gparts * sizeof(struct gpart)); */
+ /* 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]); */
+ /* } */
+ /* fclose(file); */
+
+ /* free(temp); */
+
+ error("done");
+
clocks_gettime(&time2);
e->wallclock_time = (float)clocks_diff(&time1, &time2);
diff --git a/src/runner.c b/src/runner.c
index 1a094297aaa1fd66c347ae040fe43eaab4ba4e18..6beb59737f5b603f8133d639dbf53ab437713ff8 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1135,7 +1135,6 @@ void *runner_main(void *data) {
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)
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 113019a322df6ebf4c2cdefa24d44785f744c5f1..554c60f65a3c3b653d217f793a0d971162adb351 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -58,6 +58,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
*
+ * @param ci First #cell.
+ * @param cj Second #cell.
+ *
+ * @todo Deal with periodicity.
*/
__attribute__((always_inline)) INLINE static int are_neighbours(
const struct cell *restrict ci, const struct cell *restrict cj) {
@@ -68,12 +72,12 @@ __attribute__((always_inline)) INLINE static int are_neighbours(
#endif
/* Maximum allowed distance */
- const float min_dist = ci->h[0];
+ 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 float center_i = ci->loc[k];
- const float center_j = cj->loc[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;
}
@@ -339,7 +343,12 @@ static void runner_dopair_grav(struct runner *r, struct cell *ci,
"itself.");
/* Are the cells direct neighbours? */
- if (!are_neighbours(ci, cj)) error("Non-neighbouring cells !");
+ if (!are_neighbours(ci, cj))
+ 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]);
#endif
diff --git a/src/tools.c b/src/tools.c
index db7285d917da2626d3f81e6633aefa91129f93f7..f298b3494da0c77dd76fdfb30940afa3db1a3d8e 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -507,7 +507,6 @@ void gravity_n2(struct gpart *gparts, const int gcount) {
struct gpart *restrict gpi = &gparts[pid];
const float mi = gpi->mass;
- /* Loop over every particle in the other cell. */
for (int pjd = pid + 1; pjd < gcount; pjd++) {
/* Get a hold of the jth part in ci. */