From c94fd6d695d6c67d8197142d7a2bd96a13c84d07 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sun, 1 Oct 2017 16:02:47 +0100
Subject: [PATCH] Do not escape the cells with no gparticles in the long-range
gravity task.
---
examples/EAGLE_12/eagle_12.yml | 2 +-
src/const.h | 2 +
src/engine.c | 97 +++++++++++++++++++++++---
src/gravity_properties.c | 2 +-
src/runner.c | 15 +++-
src/runner_doiact_grav.h | 121 ++++++++++++++++++++++++++-------
6 files changed, 201 insertions(+), 38 deletions(-)
diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml
index cc5d87ca02..ba181c8374 100644
--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -27,7 +27,7 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
epsilon: 0.001 # Softening length (in internal units).
- theta: 0.7 # Opening angle (Multipole acceptance criterion)
+ theta: 2.5 # Opening angle (Multipole acceptance criterion)
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/src/const.h b/src/const.h
index c8060a2be5..60c84a240a 100644
--- a/src/const.h
+++ b/src/const.h
@@ -113,4 +113,6 @@
#define SOURCETERMS_NONE
//#define SOURCETERMS_SN_FEEDBACK
+#define ICHECK 5726454604296ll
+
#endif /* SWIFT_CONST_H */
diff --git a/src/engine.c b/src/engine.c
index a773f6f4fd..b46bd79cf5 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1317,6 +1317,7 @@ void engine_addtasks_recv_gravity(struct engine *e, struct cell *c, struct task
}
c->recv_grav = t_grav;
+ c->recv_multipole = t_multi;
for (struct link *l = c->grav; l != NULL; l = l->next) {
scheduler_addunlock(s, t_grav, l->t);
@@ -1948,13 +1949,32 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
struct gravity_tensors *const multi_i = ci->multipole;
const double CoM_i[3] = {multi_i->CoM[0], multi_i->CoM[1], multi_i->CoM[2]};
- /* Loop over every other cell */
- for (int ii = 0; ii < cdim[0]; ii++) {
- for (int jj = 0; jj < cdim[1]; jj++) {
- for (int kk = 0; kk < cdim[2]; kk++) {
+ /* /\* Loop over every other cell *\/ */
+ /* for (int ii = 0; ii < cdim[0]; ii++) { */
+ /* for (int jj = 0; jj < cdim[1]; jj++) { */
+ /* for (int kk = 0; kk < cdim[2]; kk++) { */
- /* Get the cell */
- const int cjd = cell_getid(cdim, ii, jj, kk);
+ /* /\* Get the cell *\/ */
+ /* const int cjd = cell_getid(cdim, ii, jj, kk); */
+ /* struct cell *cj = &cells[cjd]; */
+
+
+ /* 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];
/* Avoid duplicates of local pairs*/
@@ -1980,9 +2000,13 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
const double r2 = dx * dx + dy * dy + dz * dz;
/* Are the cells too close for a MM interaction ? */
- if (!gravity_M2L_accept(multi_i->r_max_rebuild,
+ if (1 || !gravity_M2L_accept(multi_i->r_max_rebuild,
multi_j->r_max_rebuild, theta_crit2, r2)) {
+ if(i==11 && j==11 && k==10)
+ message("Found direct neighbour: (i,j,k)=(%d,%d,%d) (iii,jjj,kkk)=(%d,%d,%d) nodeID=%d", i,j,k, iii,jjj,kkk, cj->nodeID);
+
+
/* Ok, we need to add a direct pair calculation */
scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0,
ci, cj);
@@ -2985,6 +3009,59 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
#endif
}
+
+
+ /* Only interested in gravity tasks as of here. */
+ if (t->subtype == task_subtype_grav) {
+
+#ifdef WITH_MPI
+ /* Activate the send/recv tasks. */
+ if (ci->nodeID != engine_rank) {
+
+ /* If the local cell is active, receive data from the foreign cell. */
+ if (cj_active) {
+ scheduler_activate(s, ci->recv_grav);
+ scheduler_activate(s, ci->recv_multipole);
+ }
+
+ /* Is the foreign cell active and will need stuff from us? */
+ if (ci_active) {
+
+ struct link *l =
+ scheduler_activate_send(s, cj->send_grav, ci->nodeID);
+
+ /* 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. */
+ cell_activate_drift_gpart(l->t->ci, s);
+
+ scheduler_activate_send(s, cj->send_multipole, ci->nodeID);
+ }
+ } else if (cj->nodeID != engine_rank) {
+
+ /* If the local cell is active, receive data from the foreign cell. */
+ if (ci_active) {
+ scheduler_activate(s, cj->recv_grav);
+ scheduler_activate(s, cj->recv_multipole);
+ }
+
+ /* Is the foreign cell active and will need stuff from us? */
+ if (cj_active) {
+
+ struct link *l =
+ scheduler_activate_send(s, ci->send_grav, cj->nodeID);
+
+
+ /* 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. */
+ cell_activate_drift_gpart(l->t->ci, s);
+
+ scheduler_activate_send(s, ci->send_multipole, cj->nodeID);
+ }
+ }
+#endif
+ }
}
/* Kick/init ? */
@@ -3066,9 +3143,9 @@ void engine_print_task_counts(struct engine *e) {
int counts[task_type_count + 1];
for (int k = 0; k <= task_type_count; k++) counts[k] = 0;
for (int k = 0; k < nr_tasks; k++) {
- /* if (tasks[k].skip) */
- /* counts[task_type_count] += 1; */
- /* else */
+ if (tasks[k].skip)
+ counts[task_type_count] += 1;
+ else
counts[(int)tasks[k].type] += 1;
}
message("Total = %d (per cell = %d)", nr_tasks,
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index 27a5de0a41..cb3860b5f6 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -52,7 +52,7 @@ void gravity_props_init(struct gravity_props *p,
/* Opening angle */
p->theta_crit = parser_get_param_double(params, "Gravity:theta");
- if (p->theta_crit >= 1.) error("Theta too large. FMM won't converge.");
+ //if (p->theta_crit >= 1.) error("Theta too large. FMM won't converge.");
p->theta_crit2 = p->theta_crit * p->theta_crit;
p->theta_crit_inv = 1. / p->theta_crit;
diff --git a/src/runner.c b/src/runner.c
index 8e99a7f6cb..5afc5eca3d 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1419,6 +1419,14 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
+
+#if (ICHECK != 0)
+ for(int i=0; i < c->gcount; ++i)
+ if(c->gparts[i].id_or_neg_offset == ICHECK)
+ message("Found gpart");
+#endif
+
+
/* Anything to do here? */
if (!cell_is_active(c, e)) return;
@@ -1476,14 +1484,15 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
/* Check that this gpart has interacted with all the other
* particles (via direct or multipoles) in the box */
- if (gp->num_interacted != (long long)e->total_nr_gparts)
+ if (gp->num_interacted != e->total_nr_gparts && gp->id_or_neg_offset == ICHECK)
error(
"g-particle (id=%lld, type=%s) did not interact "
"gravitationally "
"with all other gparts gp->num_interacted=%lld, "
- "total_gparts=%zd",
+ "total_gparts=%zd (local num_gparts=%zd)",
gp->id_or_neg_offset, part_type_names[gp->type],
- gp->num_interacted, e->total_nr_gparts);
+ gp->num_interacted, e->total_nr_gparts,
+ e->s->nr_gparts);
}
#endif
}
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 8ff90412cc..89d7231d8a 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -43,6 +43,13 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
struct gpart *gparts = c->gparts;
const int gcount = c->gcount;
+#if (ICHECK != 0)
+ for(int i=0; i < c->gcount; ++i)
+ if(c->gparts[i].id_or_neg_offset == ICHECK)
+ message("Found gpart depth=%d split=%d m->num_interacted=%lld",
+ c->depth, c->split, c->multipole->pot.num_interacted);
+#endif
+
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
@@ -1136,20 +1143,40 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
const struct gravity_props *props = e->gravity_properties;
const int periodic = s->periodic;
const double cell_width = s->width[0];
+ const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const double theta_crit2 = props->theta_crit2;
const double max_distance = props->a_smooth * props->r_cut_max * cell_width;
const double max_distance2 = max_distance * max_distance;
+#if (ICHECK != 0)
+ int check = 0;
+ int direct_ngbs = 0;
+ int direct_ngbs_gpart = 0;
+ int other_ngbs_gpart = 0;
+ for(int i=0; i < ci->gcount; ++i)
+ if(ci->gparts[i].id_or_neg_offset == ICHECK) {
+ message("Found gpart");
+ check = 1;
+ }
+#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+ long long counter = 0;
+#endif
+
TIMER_TIC;
/* Recover the list of top-level cells */
- struct cell *cells = e->s->cells_top;
- const int nr_cells = e->s->nr_cells;
+ struct cell *cells = s->cells_top;
+ const int nr_cells = s->nr_cells;
/* Anything to do here? */
if (!cell_is_active(ci, e)) return;
+ if (ci->nodeID != engine_rank)
+ error("Non-local cell in long-range gravity task!");
+
/* Check multipole has been drifted */
if (ci->ti_old_multipole != e->ti_current)
error("Interacting un-drifted multipole");
@@ -1157,38 +1184,74 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
/* Recover the local multipole */
struct gravity_tensors *const multi_i = ci->multipole;
const double CoM_i[3] = {multi_i->CoM[0], multi_i->CoM[1], multi_i->CoM[2]};
- const double CoM_rebuild_i[3] = {multi_i->CoM_rebuild[0],
- multi_i->CoM_rebuild[1],
- multi_i->CoM_rebuild[2]};
+ /* const double CoM_rebuild_i[3] = {multi_i->CoM_rebuild[0], */
+ /* multi_i->CoM_rebuild[1], */
+ /* multi_i->CoM_rebuild[2]}; */
+
+ /* Get the cell index. MATTHIEU */
+ const int cid = (ci - cells);// / sizeof(struct cell);
+ const int i = cid / (cdim[1] * cdim[2]);
+ const int j = (cid / cdim[2]) % cdim[1];
+ const int k = cid % cdim[2];
/* Loop over all the top-level cells and go for a M-M interaction if
* well-separated */
- for (int i = 0; i < nr_cells; ++i) {
+ for (int n = 0; n < nr_cells; ++n) {
/* Handle on the top-level cell and it's gravity business*/
- struct cell *cj = &cells[i];
+ struct cell *cj = &cells[n];
const struct gravity_tensors *const multi_j = cj->multipole;
- /* Avoid stupid cases */
- if (ci == cj || cj->gcount == 0) continue;
+ /* Avoid self contributions */
+ if (ci == cj) continue;
- /* Get the distance between the CoMs at the last rebuild*/
- double dx_r = CoM_rebuild_i[0] - multi_j->CoM_rebuild[0];
- double dy_r = CoM_rebuild_i[1] - multi_j->CoM_rebuild[1];
- double dz_r = CoM_rebuild_i[2] - multi_j->CoM_rebuild[2];
+#ifdef SWIFT_DEBUG_CHECKS
+ counter += multi_j->m_pole.num_gpart;
+#endif
- /* Apply BC */
- if (periodic) {
- dx_r = nearest(dx_r, dim[0]);
- dy_r = nearest(dy_r, dim[1]);
- dz_r = nearest(dz_r, dim[2]);
- }
- const double r2_rebuild = dx_r * dx_r + dy_r * dy_r + dz_r * dz_r;
+ // MATTHIEU
+ const int cjd = (cj - cells);// / sizeof(struct cell);
+ const int ii = cjd / (cdim[1] * cdim[2]);
+ const int jj = (cjd / cdim[2]) % cdim[1];
+ const int kk = cjd % cdim[2];
+ /* if(check) message("cid=%d cjd=%d cj->nodeID=%d cj->gcount=%d", */
+ /* cid, cjd, cj->nodeID, cj->gcount); */
+
+ /* /\* Get the distance between the CoMs at the last rebuild*\/ */
+ /* double dx_r = CoM_rebuild_i[0] - multi_j->CoM_rebuild[0]; */
+ /* double dy_r = CoM_rebuild_i[1] - multi_j->CoM_rebuild[1]; */
+ /* double dz_r = CoM_rebuild_i[2] - multi_j->CoM_rebuild[2]; */
+
+ /* /\* Apply BC *\/ */
+ /* if (periodic) { */
+ /* dx_r = nearest(dx_r, dim[0]); */
+ /* dy_r = nearest(dy_r, dim[1]); */
+ /* dz_r = nearest(dz_r, dim[2]); */
+ /* } */
+ /* const double r2_rebuild = dx_r * dx_r + dy_r * dy_r + dz_r * dz_r; */
+
+ /* Are we in charge of this cell pair? MATTHIEU*/
+ /* if (gravity_M2L_accept(multi_i->r_max_rebuild, multi_j->r_max_rebuild, */
+ /* theta_crit2, r2_rebuild)) { */
+ if((abs(i-ii) <= 1 || abs(i-ii - cdim[0]) <= 1) &&
+ (abs(j-jj) <= 1 || abs(j-jj - cdim[1]) <= 1) &&
+ (abs(k-kk) <= 1)|| abs(k-kk - cdim[2]) <= 1) {
+
+#if (ICHECK != 0)
+ if(check) {
+ ++direct_ngbs;
+ direct_ngbs_gpart += cj->multipole->m_pole.num_gpart;
+ message("Found direct neighbour %d: (i,j,k)=(%d,%d,%d) (ii,jj,kk)=(%d,%d,%d) nodeID=%d",
+ direct_ngbs, i,j,k, ii,jj,kk, cj->nodeID);
+ }
+#endif
+
- /* Are we in charge of this cell pair? */
- if (gravity_M2L_accept(multi_i->r_max_rebuild, multi_j->r_max_rebuild,
- theta_crit2, r2_rebuild)) {
+ }else{
+ if(check)
+ other_ngbs_gpart += cj->multipole->m_pole.num_gpart;
+
/* Let's compute the current distance between the cell pair*/
double dx = CoM_i[0] - multi_j->CoM[0];
double dy = CoM_i[1] - multi_j->CoM[1];
@@ -1225,6 +1288,18 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
} /* We are in charge of this pair */
} /* Loop over top-level cells */
+
+ if(check)
+ message("Interacted with %d indirectly and ignored %d direct interactions (counter=%lld) nr_cells=%d",
+ other_ngbs_gpart, direct_ngbs_gpart, counter, nr_cells);
+
+
+#ifdef SWIFT_DEBUG_CHECKS
+ counter += ci->multipole->m_pole.num_gpart;
+ if(counter != e->total_nr_gparts)
+ error("Not found the right number of particles in top-level interactions");
+#endif
+
if (timer) TIMER_TOC(timer_dograv_long_range);
}
--
GitLab