diff --git a/src/active.h b/src/active.h
index d53c2fcce750e4ea2d9b280cc968f0dc1dd04304..fedb668c8f14bb7c3dade20990a85fb39f3bf68b 100644
--- a/src/active.h
+++ b/src/active.h
@@ -129,7 +129,7 @@ __attribute__((always_inline)) INLINE static int cell_is_active_gravity(
const struct cell *c, const struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
- if (c->ti_gravity_end_min < e->ti_current && c->nr_tasks > 0)
+ if (c->ti_gravity_end_min < e->ti_current)
error(
"cell in an impossible time-zone! c->ti_end_min=%lld (t=%e) and "
"e->ti_current=%lld (t=%e, a=%e)",
@@ -140,6 +140,20 @@ __attribute__((always_inline)) INLINE static int cell_is_active_gravity(
return (c->ti_gravity_end_min == e->ti_current);
}
+/**
+ * @brief Does a cell contain any g-particle finishing their time-step now ?
+ *
+ * @param c The #cell.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #cell contains at least an active particle, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int cell_is_active_gravity_mm(
+ const struct cell *c, const struct engine *e) {
+
+ return (c->ti_gravity_end_min == e->ti_current);
+}
+
+
/**
* @brief Are *all* g-particles in a cell finishing their time-step now ?
*
diff --git a/src/cell.c b/src/cell.c
index 1ceb0556a5bfada87d2a5b160dd6c353724dc315..4d7cf05b92d4020e3d56fccbe9fe53d2fe0e23c3 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -67,6 +67,10 @@
/* Global variables. */
int cell_next_tag = 0;
+integertime_t ti_current_global;
+
+extern FILE* file_dump;
+
/**
* @brief Get the size of the cell subtree.
*
@@ -180,9 +184,12 @@ int cell_pack(struct cell *restrict c, struct pcell *restrict pc) {
pc->ti_hydro_end_max = c->ti_hydro_end_max;
pc->ti_gravity_end_min = c->ti_gravity_end_min;
pc->ti_gravity_end_max = c->ti_gravity_end_max;
- pc->ti_old_part = c->ti_old_part;
- pc->ti_old_gpart = c->ti_old_gpart;
- pc->ti_old_multipole = c->ti_old_multipole;
+ //pc->ti_old_part = c->ti_old_part;
+ // pc->ti_old_gpart = c->ti_old_gpart;
+ //pc->ti_old_multipole = c->ti_old_multipole;
+
+ pc->multipole = *(c->multipole);
+
pc->count = c->count;
pc->gcount = c->gcount;
pc->scount = c->scount;
@@ -230,9 +237,16 @@ int cell_unpack(struct pcell *restrict pc, struct cell *restrict c,
c->ti_hydro_end_max = pc->ti_hydro_end_max;
c->ti_gravity_end_min = pc->ti_gravity_end_min;
c->ti_gravity_end_max = pc->ti_gravity_end_max;
- c->ti_old_part = pc->ti_old_part;
- c->ti_old_gpart = pc->ti_old_gpart;
- c->ti_old_multipole = pc->ti_old_multipole;
+ //c->ti_old_part = pc->ti_old_part;
+ //c->ti_old_gpart = pc->ti_old_gpart;
+ //c->ti_old_multipole = pc->ti_old_multipole;
+
+ c->ti_old_part = ti_current_global;
+ c->ti_old_gpart = ti_current_global;
+ c->ti_old_multipole = ti_current_global;
+
+ *(c->multipole) = pc->multipole;
+
c->count = pc->count;
c->gcount = pc->gcount;
c->scount = pc->scount;
@@ -339,6 +353,11 @@ int cell_unpack_end_step(struct cell *restrict c,
c->ti_gravity_end_max = pcells[0].ti_gravity_end_max;
c->dx_max_part = pcells[0].dx_max_part;
+ if(c->ti_last_update == ti_current_global && c->ti_last_update > 0)
+ error("Overwriting data");
+
+ c->ti_last_update = ti_current_global;
+
/* Fill in the progeny, depth-first recursion. */
int count = 1;
for (int k = 0; k < 8; k++)
@@ -1878,7 +1897,7 @@ void cell_activate_grav_mm_task(struct cell *ci, struct cell *cj,
const struct engine *e = s->space->e;
/* Anything to do here? */
- if (!cell_is_active_gravity(ci, e) && !cell_is_active_gravity(cj, e))
+ if (!cell_is_active_gravity_mm(ci, e) && !cell_is_active_gravity_mm(cj, e))
error("Inactive MM task being activated");
/* Atomically drift the multipole in ci */
@@ -2273,7 +2292,9 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
} else if (t->type == task_type_pair) {
cell_activate_subcell_grav_tasks(ci, cj, s);
} else if (t->type == task_type_grav_mm) {
- error("Incorrectyl linked M-M task!");
+#ifdef SWIFT_DEBUG_CHECKS
+ error("Incorrectly linked M-M task!");
+#endif
}
}
@@ -2284,8 +2305,18 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
if (ci_nodeID != nodeID) {
/* If the local cell is active, receive data from the foreign cell. */
- if (cj_active) scheduler_activate(s, ci->recv_grav);
+ if (cj_active) {scheduler_activate(s, ci->recv_grav);
+
+ fprintf(file_dump, "node %d activating cell %d to be received. ",
+ cj->nodeID, ci->cellID);
+ if(ci->parent == NULL)
+ fprintf(file_dump, "top\n");
+ else
+ fprintf(file_dump, " \n");
+ fflush(file_dump);
+
+ }
/* If the foreign cell is active, we want its ti_end values. */
if (ci_active) scheduler_activate(s, ci->recv_ti);
@@ -2294,6 +2325,15 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate_send(s, cj->send_grav, ci_nodeID);
+ fprintf(file_dump, "node %d activating cell %d to send to node %d ",
+ cj->nodeID, cj->cellID, ci->nodeID);
+ if(cj->parent == NULL)
+ fprintf(file_dump, "top\n");
+ else
+ fprintf(file_dump, " \n");
+
+ fflush(file_dump);
+
/* Drift the cell which will be sent at the level at which it is
sent, i.e. drift the cell specified in the send task (l->t)
itself. */
@@ -2306,7 +2346,18 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
} else if (cj_nodeID != nodeID) {
/* If the local cell is active, receive data from the foreign cell. */
- if (ci_active) scheduler_activate(s, cj->recv_grav);
+ if (ci_active) {scheduler_activate(s, cj->recv_grav);
+
+ fprintf(file_dump, "node %d activating cell %d to be received. ",
+ ci->nodeID, cj->cellID);
+ if(cj->parent == NULL)
+ fprintf(file_dump, "top\n");
+ else
+ fprintf(file_dump, " \n");
+
+
+ fflush(file_dump);
+ }
/* If the foreign cell is active, we want its ti_end values. */
if (cj_active) scheduler_activate(s, cj->recv_ti);
@@ -2316,6 +2367,15 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate_send(s, ci->send_grav, cj_nodeID);
+ fprintf(file_dump, "node %d activating cell %d to sent to node %d",
+ ci->nodeID, ci->cellID, cj->nodeID);
+ if(ci->parent == NULL)
+ fprintf(file_dump, " top\n");
+ else
+ fprintf(file_dump, " \n");
+
+ fflush(file_dump);
+
/* Drift the cell which will be sent at the level at which it is
sent, i.e. drift the cell specified in the send task (l->t)
itself. */
@@ -2334,8 +2394,8 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
struct task *t = l->t;
struct cell *ci = t->ci;
struct cell *cj = t->cj;
- const int ci_active = cell_is_active_gravity(ci, e);
- const int cj_active = (cj != NULL) ? cell_is_active_gravity(cj, e) : 0;
+ const int ci_active = cell_is_active_gravity_mm(ci, e);
+ const int cj_active = (cj != NULL) ? cell_is_active_gravity_mm(cj, e) : 0;
#ifdef WITH_MPI
const int ci_nodeID = ci->nodeID;
const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
@@ -2344,7 +2404,9 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
const int cj_nodeID = nodeID;
#endif
+#ifdef SWIFT_DEBUG_CHECKS
if (t->type != task_type_grav_mm) error("Incorrectly linked gravity task!");
+#endif
/* Only activate tasks that involve a local active cell. */
if ((ci_active && ci_nodeID == nodeID) ||
@@ -2356,7 +2418,7 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
}
/* Unskip all the other task types. */
- if (c->nodeID == nodeID && cell_is_active_gravity(c, e)) {
+ if (c->nodeID == nodeID && cell_is_active_gravity_mm(c, e)) {
if (c->init_grav != NULL) scheduler_activate(s, c->init_grav);
if (c->init_grav_out != NULL) scheduler_activate(s, c->init_grav_out);
@@ -2898,3 +2960,57 @@ int cell_can_use_pair_mm(const struct cell *ci, const struct cell *cj,
return gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2);
}
+
+
+/**
+ * @brief Can we use the MM interactions fo a given pair of cells?
+ *
+ * @param ci The first #cell.
+ * @param cj The second #cell.
+ * @param e The #engine.
+ * @param s The #space.
+ */
+int cell_can_use_pair_mm_rebuild(const struct cell *ci, const struct cell *cj,
+ const struct engine *e, const struct space *s) {
+
+ const double theta_crit2 = e->gravity_properties->theta_crit2;
+ const int periodic = s->periodic;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
+
+ /* Recover the multipole information */
+ const struct gravity_tensors *const multi_i = ci->multipole;
+ const struct gravity_tensors *const multi_j = cj->multipole;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ if(multi_i->CoM[0] < ci->loc[0] || multi_i->CoM[0] > ci->loc[0] + ci->width[0])
+ error("Invalid multipole position ci");
+ if(multi_i->CoM[1] < ci->loc[1] || multi_i->CoM[1] > ci->loc[1] + ci->width[1])
+ error("Invalid multipole position ci");
+ if(multi_i->CoM[2] < ci->loc[2] || multi_i->CoM[2] > ci->loc[2] + ci->width[2])
+ error("Invalid multipole position ci");
+
+ if(multi_j->CoM[0] < cj->loc[0] || multi_j->CoM[0] > cj->loc[0] + cj->width[0])
+ error("Invalid multipole position cj");
+ if(multi_j->CoM[1] < cj->loc[1] || multi_j->CoM[1] > cj->loc[1] + cj->width[1])
+ error("Invalid multipole position cj");
+ if(multi_j->CoM[2] < cj->loc[2] || multi_j->CoM[2] > cj->loc[2] + cj->width[2])
+ error("Invalid multipole position cj");
+
+#endif
+
+ /* Get the distance between the CoMs */
+ double dx = multi_i->CoM[0] - multi_j->CoM[0];
+ double dy = multi_i->CoM[1] - multi_j->CoM[1];
+ double dz = multi_i->CoM[2] - multi_j->CoM[2];
+
+ /* Apply BC */
+ if (periodic) {
+ dx = nearest(dx, dim[0]);
+ dy = nearest(dy, dim[1]);
+ dz = nearest(dz, dim[2]);
+ }
+ const double r2 = dx * dx + dy * dy + dz * dz;
+
+ return gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2);
+}
diff --git a/src/cell.h b/src/cell.h
index ed017d2d11f6612321df5a7bb9c514dd3b5179f9..fc18ed510ab20ea754d6be887f05a86a20fe8cf5 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -122,6 +122,8 @@ struct pcell {
/*! Relative indices of the cell's progeny. */
int progeny[8];
+ struct gravity_tensors multipole;
+
#ifdef SWIFT_DEBUG_CHECKS
/* Cell ID (for debugging) */
int cellID;
@@ -157,6 +159,8 @@ struct pcell_step {
*/
struct cell {
+ integertime_t ti_last_update;
+
/*! The cell location on the grid. */
double loc[3];
@@ -529,6 +533,8 @@ void cell_set_super_mapper(void *map_data, int num_elements, void *extra_data);
int cell_has_tasks(struct cell *c);
int cell_can_use_pair_mm(const struct cell *ci, const struct cell *cj,
const struct engine *e, const struct space *s);
+int cell_can_use_pair_mm_rebuild(const struct cell *ci, const struct cell *cj,
+ const struct engine *e, const struct space *s);
/* Inlined functions (for speed). */
diff --git a/src/engine.c b/src/engine.c
index 78705411d455ab032d184410cad3f957b7b0cfec..f35e4d80b16bae013e2737e87c13ab428cafffc2 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -92,6 +92,8 @@
/* Particle cache size. */
#define CACHE_SIZE 512
+FILE* file_dump;
+
const char *engine_policy_names[] = {"none",
"rand",
"steal",
@@ -114,6 +116,8 @@ const char *engine_policy_names[] = {"none",
/** The rank of the engine as a global variable (for messages). */
int engine_rank;
+extern integertime_t ti_current_global;
+
/**
* @brief Data collected from the cells at the end of a time-step
*/
@@ -2496,7 +2500,7 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
if (periodic && min_radius > max_distance) continue;
/* Are the cells too close for a MM interaction ? */
- if (!cell_can_use_pair_mm(ci, cj, e, s)) {
+ if (!cell_can_use_pair_mm_rebuild(ci, cj, e, s)) {
/* Ok, we need to add a direct pair calculation */
scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0,
@@ -3739,8 +3743,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
const struct cell *cj = t->cj;
const int ci_nodeID = ci->nodeID;
const int cj_nodeID = cj->nodeID;
- const int ci_active_gravity = cell_is_active_gravity(ci, e);
- const int cj_active_gravity = cell_is_active_gravity(cj, e);
+ const int ci_active_gravity = cell_is_active_gravity_mm(ci, e);
+ const int cj_active_gravity = cell_is_active_gravity_mm(cj, e);
if ((ci_active_gravity && ci_nodeID == engine_rank) ||
(cj_active_gravity && cj_nodeID == engine_rank))
@@ -4003,12 +4007,12 @@ void engine_rebuild(struct engine *e, int clean_smoothing_length_values) {
/* If in parallel, exchange the cell structure, top-level and neighbouring
* multipoles. */
#ifdef WITH_MPI
- engine_exchange_cells(e);
-
if (e->policy & engine_policy_self_gravity) engine_exchange_top_multipoles(e);
- if (e->policy & engine_policy_self_gravity)
- engine_exchange_proxy_multipoles(e);
+ engine_exchange_cells(e);
+
+ //if (e->policy & engine_policy_self_gravity)
+ // engine_exchange_proxy_multipoles(e);
#endif
/* Re-build the tasks. */
@@ -4761,7 +4765,11 @@ void engine_step(struct engine *e) {
e->min_active_bin = get_min_active_bin(e->ti_current, e->ti_old);
e->step += 1;
e->step_props = engine_step_prop_none;
+
+ ti_current_global = e->ti_current;
+ fprintf(file_dump, "#################### Step %d #####################\n\n\n", e->step);
+
/* When restarting, move everyone to the current time. */
if (e->restarting) engine_drift_all(e);
@@ -5344,7 +5352,7 @@ void engine_makeproxies(struct engine *e) {
if (e->verbose)
message(
"Looking for proxies up to %d top-level cells away (delta_m=%d "
- "delta_m=%d)",
+ "delta_p=%d)",
delta, delta_m, delta_p);
/* Loop over each cell in the space. */
@@ -5413,7 +5421,16 @@ void engine_makeproxies(struct engine *e) {
if (with_gravity) {
/* Are we too close for M2L? */
- if (!cell_can_use_pair_mm(&cells[cid], &cells[cjd], e, s))
+ //if (!cell_can_use_pair_mm_rebuild(&cells[cid], &cells[cjd], e, s))
+ /* if (((abs(ind[0] - ii) <= 5 || */
+ /* abs(ind[0] - ii - cdim[0]) <= 5 || */
+ /* abs(ind[0] - ii + cdim[0]) <= 5) && */
+ /* (abs(ind[1] - jj) <= 5 || */
+ /* abs(ind[1] - jj - cdim[1]) <= 5 || */
+ /* abs(ind[1] - jj + cdim[1]) <= 5) && */
+ /* (abs(ind[2] - kk) <= 5 || */
+ /* abs(ind[2] - kk - cdim[2]) <= 5 || */
+ /* abs(ind[2] - kk + cdim[2]) <= 5))) */
proxy_type |= (int)proxy_cell_type_gravity;
}
@@ -5893,6 +5910,10 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
e->timeFirstSTFOutput = 0;
engine_rank = nodeID;
+ char buffer[32];
+ sprintf(buffer, "dump_%d.txt", engine_rank);
+ file_dump = fopen(buffer, "w");
+
/* Initialise VELOCIraptor. */
if (e->policy & engine_policy_structure_finding) {
parser_get_param_string(params, "StructureFinding:basename",
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 812a6e3850fa66f614326c872511c36ba91b5635..69213095add9b9f4bf4425ddf0cd3f723497deb3 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -1144,7 +1144,7 @@ static INLINE void runner_dopair_grav_mm(struct runner *r,
TIMER_TIC;
/* Anything to do here? */
- if (!cell_is_active_gravity(ci, e) || ci->nodeID != engine_rank) return;
+ if (!cell_is_active_gravity_mm(ci, e) || ci->nodeID != engine_rank) return;
/* Short-cut to the multipole */
const struct multipole *multi_j = &cj->multipole->m_pole;
@@ -1485,8 +1485,8 @@ static INLINE void runner_dopair_grav_mm_symmetric(struct runner *r,
error("Running M-M task with two inactive cells.");
#endif
- if (cell_is_active_gravity(ci, e)) runner_dopair_grav_mm(r, ci, cj);
- if (cell_is_active_gravity(cj, e)) runner_dopair_grav_mm(r, cj, ci);
+ if (cell_is_active_gravity_mm(ci, e)) runner_dopair_grav_mm(r, ci, cj);
+ if (cell_is_active_gravity_mm(cj, e)) runner_dopair_grav_mm(r, cj, ci);
}
/**
diff --git a/src/scheduler.c b/src/scheduler.c
index 16e627f27a3d186506105a56249f9d3f4a8fe8aa..78bbdcfb042350afcb0e5d65739f2258351027ba 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -889,7 +889,7 @@ static void scheduler_splittask_gravity(struct task *t, struct scheduler *s) {
}
/* Should we replace it with an M-M task? */
- if (cell_can_use_pair_mm(ci, cj, e, sp)) {
+ if (cell_can_use_pair_mm_rebuild(ci, cj, e, sp)) {
t->type = task_type_grav_mm;
t->subtype = task_subtype_none;