Commit ebb53467 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Record r_max and CoM of each multipole at rebuild time.

parent 107aa7e5
......@@ -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 */
......@@ -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(
......
......@@ -638,7 +638,7 @@ 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,
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);
......
......@@ -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. */
......
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