diff --git a/src/multipole.h b/src/multipole.h index 48be0a31dd2d389f0fa742feb15eec5fc0b11646..15f9371fa2e8eda14b94cac031f153b498d066a8 100644 --- a/src/multipole.h +++ b/src/multipole.h @@ -176,9 +176,15 @@ struct gravity_tensors { /*! Centre of mass of the matter dsitribution */ double CoM[3]; + /*! Centre of mass of the matter dsitribution at the last rebuild */ + double CoM_rebuild[3]; + /*! Upper limit of the CoM<->gpart distance */ double r_max; + /*! Upper limit of the CoM<->gpart distance at the last rebuild */ + double r_max_rebuild; + /*! Multipole mass */ struct multipole m_pole; @@ -2537,4 +2543,30 @@ __attribute__((always_inline)) INLINE static int gravity_multipole_accept( return (r2 > (r_crit_a + r_crit_b) * (r_crit_a + r_crit_b)); } +/** + * @brief Checks whether a cell-cell interaction can be appromixated by a M-M + * interaction using the r_max at rebuild time. + * + * @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_rebuild(const struct gravity_tensors *ma, + const struct gravity_tensors *mb, + double theta_crit_inv) { + + const double r_crit_a = ma->r_max_rebuild * theta_crit_inv; + const double r_crit_b = mb->r_max_rebuild * theta_crit_inv; + + const double dx = ma->CoM_rebuild[0] - mb->CoM_rebuild[0]; + const double dy = ma->CoM_rebuild[1] - mb->CoM_rebuild[1]; + const double dz = ma->CoM_rebuild[2] - mb->CoM_rebuild[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.c b/src/runner.c index 18e48c945db389e935924b3356d68d63349f45cf..512a4f362fe80caf0441c191f5503e94aec7e763 100644 --- a/src/runner.c +++ b/src/runner.c @@ -1392,8 +1392,7 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) { if (e->policy & engine_policy_self_gravity) { /* Check that this gpart has interacted with all the other - * particles - * (via direct or multipoles) in the box */ + * particles (via direct or multipoles) in the box */ gp->num_interacted++; if (gp->num_interacted != (long long)e->s->nr_gparts) error( diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h index 8161a13e765b032642f63c81a0c3c9b27e35d51b..d35e100ec63a1a647ec1e145e65bf0c672c19d35 100644 --- a/src/runner_doiact_grav.h +++ b/src/runner_doiact_grav.h @@ -638,8 +638,8 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) { if (ci == cj || cj->gcount == 0) continue; /* Check the multipole acceptance criterion */ - if (gravity_multipole_accept(ci->multipole, cj->multipole, - theta_crit_inv)) { + if (gravity_multipole_accept_rebuild(ci->multipole, cj->multipole, + theta_crit_inv)) { /* Go for a (non-symmetric) M-M calculation */ runner_dopair_grav_mm(r, ci, cj); } diff --git a/src/space.c b/src/space.c index 7b59be248d7a839500c70c237f22734f60264ffc..8e5c6ba15a9e4a3416f508519c85d04ca8e2a83d 100644 --- a/src/space.c +++ b/src/space.c @@ -2115,6 +2115,10 @@ void space_split_recursive(struct space *s, struct cell *c, /* Take minimum of both limits */ c->multipole->r_max = min(r_max, sqrt(dx * dx + dy * dy + dz * dz)); + c->multipole->r_max_rebuild = c->multipole->r_max; + c->multipole->CoM_rebuild[0] = c->multipole->CoM[0]; + c->multipole->CoM_rebuild[1] = c->multipole->CoM[1]; + c->multipole->CoM_rebuild[2] = c->multipole->CoM[2]; } /* Deal with gravity */ } @@ -2194,6 +2198,10 @@ void space_split_recursive(struct space *s, struct cell *c, c->multipole->r_max = 0.; } } + c->multipole->r_max_rebuild = c->multipole->r_max; + c->multipole->CoM_rebuild[0] = c->multipole->CoM[0]; + c->multipole->CoM_rebuild[1] = c->multipole->CoM[1]; + c->multipole->CoM_rebuild[2] = c->multipole->CoM[2]; } /* Set the values for this cell. */