Skip to content
GitLab
Commit c94fd6d6 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Do not escape the cells with no gparticles in the long-range gravity task.

parent cf7019f8
Hide whitespace changes
Inline Side-by-side
examples/EAGLE_12/eagle_12.yml
View file @ c94fd6d6
......@@ -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:
......
src/const.h
View file @ c94fd6d6
......@@ -113,4 +113,6 @@
#define SOURCETERMS_NONE
//#define SOURCETERMS_SN_FEEDBACK
#define ICHECK 5726454604296ll
#endif /* SWIFT_CONST_H */
src/engine.c
View file @ c94fd6d6
......@@ -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,
......
src/gravity_properties.c
View file @ c94fd6d6
......@@ -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;
......
src/runner.c
View file @ c94fd6d6
......@@ -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
}
......
src/runner_doiact_grav.h
View file @ c94fd6d6
......@@ -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);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment