diff --git a/src/engine.c b/src/engine.c index c8ce267982adbf67349fb9fc81fe08989b0b12a4..7dc6c854aacf16262044381cdcf44844550aa6b9 100644 --- a/src/engine.c +++ b/src/engine.c @@ -1647,9 +1647,9 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts, * @brief Constructs the top-level tasks for the short-range gravity * interactions. * - * All top-cells get a self task. - * All neighbouring pairs get a pair task. - * All non-neighbouring pairs within a range of 6 cells get a M-M task. + * - All top-cells get a self task. + * - All pairs within range according to the multipole acceptance + * criterion get a pair task. * * @param e The #engine. */ @@ -1658,6 +1658,7 @@ void engine_make_self_gravity_tasks(struct engine *e) { struct space *s = e->s; struct scheduler *sched = &e->sched; const int nodeID = e->nodeID; + const double theta_crit_inv = e->gravity_properties->theta_crit_inv; struct cell *cells = s->cells_top; const int nr_cells = s->nr_cells; @@ -1668,13 +1669,14 @@ void engine_make_self_gravity_tasks(struct engine *e) { /* Skip cells without gravity particles */ if (ci->gcount == 0) continue; - /* Is that neighbour local ? */ + /* Is that cell local ? */ if (ci->nodeID != nodeID) continue; /* If the cells is local build a self-interaction */ scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, NULL, 0); + /* Loop over every other cell */ for (int cjd = cid + 1; cjd < nr_cells; ++cjd) { struct cell *cj = &cells[cjd]; @@ -1683,9 +1685,11 @@ void engine_make_self_gravity_tasks(struct engine *e) { if (cj->gcount == 0) continue; /* Is that neighbour local ? */ - if (cj->nodeID != nodeID) continue; + if (cj->nodeID != nodeID) continue; // MATTHIEU - if (cell_are_neighbours(ci, cj)) + /* Are the cells to close for a MM interaction ? */ + if (!gravity_multipole_accept(ci->multipole, cj->multipole, + theta_crit_inv)) scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci, cj, 1); } diff --git a/src/gravity_properties.c b/src/gravity_properties.c index 03825c9f714de2d6aa67ca60f52970bb8e345666..f52029fa1543ad1f8d0121c8c4e6d362227f4c53 100644 --- a/src/gravity_properties.c +++ b/src/gravity_properties.c @@ -48,8 +48,8 @@ void gravity_props_init(struct gravity_props *p, p->eta = parser_get_param_float(params, "Gravity:eta"); /* Opening angle */ - p->theta = parser_get_param_double(params, "Gravity:theta"); - p->theta_inv = 1. / p->theta; + p->theta_crit = parser_get_param_double(params, "Gravity:theta"); + p->theta_crit_inv = 1. / p->theta_crit; /* Softening lengths */ p->epsilon = parser_get_param_double(params, "Gravity:epsilon"); @@ -63,7 +63,7 @@ void gravity_props_print(const struct gravity_props *p) { message("Self-gravity time integration: eta=%.4f", p->eta); - message("Self-gravity opening angle: theta=%.4f", p->theta); + message("Self-gravity opening angle: theta=%.4f", p->theta_crit); message("Self-gravity softening: epsilon=%.4f", p->epsilon); @@ -80,7 +80,7 @@ void gravity_props_print_snapshot(hid_t h_grpgrav, io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta); io_write_attribute_f(h_grpgrav, "Softening length", p->epsilon); - io_write_attribute_f(h_grpgrav, "Opening angle", p->theta); + io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit); io_write_attribute_f(h_grpgrav, "MM a_smooth", p->a_smooth); io_write_attribute_f(h_grpgrav, "MM r_cut", p->r_cut); } diff --git a/src/gravity_properties.h b/src/gravity_properties.h index a8402aab50dc7e2da8895c12bb767908e1cc9a8b..be26f0d1d23b8cec71fa3cbbeedac9f61f337b2c 100644 --- a/src/gravity_properties.h +++ b/src/gravity_properties.h @@ -42,10 +42,10 @@ struct gravity_props { float eta; /*! Tree opening angle (Multipole acceptance criterion) */ - double theta; + double theta_crit; /*! Inverse of opening angle */ - double theta_inv; + double theta_crit_inv; /*! Softening length */ double epsilon; diff --git a/src/multipole.h b/src/multipole.h index c1d9d605e08aa044750c7e515d268c7f47c0631e..e7864c97de1956c77c1f07010a868fe71f3a3e4e 100644 --- a/src/multipole.h +++ b/src/multipole.h @@ -208,9 +208,9 @@ INLINE static void gravity_reset(struct gravity_tensors *m) { INLINE static void gravity_drift(struct gravity_tensors *m, double dt) { /* Move the whole thing according to bulk motion */ - m->CoM[0] += m->m_pole.vel[0] * dt; - m->CoM[1] += m->m_pole.vel[1] * dt; - m->CoM[2] += m->m_pole.vel[2] * dt; + m->CoM[0] += m->m_pole.vel[0] * dt * 0; // MATTHIEU + m->CoM[1] += m->m_pole.vel[1] * dt * 0; + m->CoM[2] += m->m_pole.vel[2] * dt * 0; } /** @@ -2508,4 +2508,29 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb, gp->a_grav[2] += a_grav[2]; } +/** + * @brief Checks whether a cell-cell interaction can be appromixated by a M-M + * interaction. + * + * @param ma The #multipole of the first #cell. + * @param mb The #multipole of the second #cell. + * @param theta_crit_inv The inverse of the critical opening angle. + */ +__attribute__((always_inline)) INLINE static int gravity_multipole_accept( + const struct gravity_tensors *ma, const struct gravity_tensors *mb, + double theta_crit_inv) { + + const double r_crit_a = ma->r_max * theta_crit_inv; + const double r_crit_b = mb->r_max * theta_crit_inv; + + const double dx = ma->CoM[0] - mb->CoM[0]; + const double dy = ma->CoM[1] - mb->CoM[1]; + const double dz = ma->CoM[2] - mb->CoM[2]; + + const double r2 = dx * dx + dy * dy + dz * dz; + + /* Multipole acceptance criterion (Dehnen 2002, eq.10) */ + return (r2 > (r_crit_a + r_crit_b) * (r_crit_a + r_crit_b)); +} + #endif /* SWIFT_MULTIPOLE_H */ diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h index 9af2dc1105172ff259caa27de9c19d69ce316dba..e67200ecba9fc315c2106372dd6135d542c7dc3e 100644 --- a/src/runner_doiact_grav.h +++ b/src/runner_doiact_grav.h @@ -23,6 +23,7 @@ /* Includes. */ #include "cell.h" #include "gravity.h" +#include "inline.h" #include "part.h" /** @@ -379,14 +380,15 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj, "The impossible has happened: pair interaction between a cell and " "itself."); - /* Are the cells direct neighbours? */ - if (!cell_are_neighbours(ci, cj)) - error( - "Non-neighbouring cells ! ci->x=[%f %f %f] ci->width=%f cj->loc=[%f %f " - "%f] " - "cj->width=%f", - ci->loc[0], ci->loc[1], ci->loc[2], ci->width[0], cj->loc[0], - cj->loc[1], cj->loc[2], cj->width[0]); +/* Are the cells direct neighbours? */ +/* if (!cell_are_neighbours(ci, cj)) */ +/* error( */ +/* "Non-neighbouring cells ! ci->x=[%f %f %f] ci->width=%f cj->loc=[%f %f + * " */ +/* "%f] " */ +/* "cj->width=%f", */ +/* ci->loc[0], ci->loc[1], ci->loc[2], ci->width[0], cj->loc[0], */ +/* cj->loc[1], cj->loc[2], cj->width[0]); */ #endif @@ -515,6 +517,14 @@ void runner_dosub_grav(struct runner *r, struct cell *ci, struct cell *cj, } } +/** + * @brief Performs all M-M interactions between a given top-level cell and all + * the other top-levels that are far enough. + * + * @param r The thread #runner. + * @param ci The #cell of interest. + * @param timer Are we timing this ? + */ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) { #if ICHECK > 0 @@ -530,39 +540,40 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) { } #endif + /* Some constants */ + const struct engine *e = r->e; + const struct gravity_props *props = e->gravity_properties; + const double theta_crit_inv = props->theta_crit_inv; + TIMER_TIC; /* Recover the list of top-level cells */ - const struct engine *e = r->e; struct cell *cells = e->s->cells_top; const int nr_cells = e->s->nr_cells; - /* const double max_d = */ - /* const_gravity_a_smooth * const_gravity_r_cut * ci->width[0]; */ - /* const double max_d2 = max_d * max_d; */ - // const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]}; /* Anything to do here? */ - if (!cell_is_active(ci, e)) return; + if (!cell_is_active(ci, e)) return; // MATTHIEU (should never happen) - /* Drift our own multipole if need be */ - if (ci->ti_old_multipole != e->ti_current) cell_drift_multipole(ci, e); + /* Check multipole has been drifted */ + if (ci->ti_old_multipole != e->ti_current) + error("Interacting un-drifted multipole"); - /* Loop over all the cells and go for a p-m interaction if far enough but not - * too far */ + /* 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) { + /* Handle on the top-level cell */ struct cell *cj = &cells[i]; - if (ci == cj) continue; - if (cj->gcount == 0) continue; - - /* const double dx[3] = {cj->loc[0] - pos_i[0], // x */ - /* cj->loc[1] - pos_i[1], // y */ - /* cj->loc[2] - pos_i[2]}; // z */ - /* const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; */ - /* if (r2 > max_d2) continue; */ + /* Avoid stupid cases */ + if (ci == cj || cj->gcount == 0) continue; - if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_mm(r, ci, cj); + /* Check the multipole acceptance criterion */ + if (gravity_multipole_accept(ci->multipole, cj->multipole, + theta_crit_inv)) { + /* Go for a (non-symmetric) M-M calculation */ + runner_dopair_grav_mm(r, ci, cj); + } } if (timer) TIMER_TOC(timer_dograv_long_range);