Commit 1b4cc45e authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Branch to a different, more accurate, calculation in the long-range that don't...

Branch to a different, more accurate, calculation in the long-range that don't match the MAC any more.
parent ebb53467
......@@ -1693,7 +1693,7 @@ void engine_make_self_gravity_tasks(struct engine *e) {
/* Are the cells to close for a MM interaction ? */
if (!gravity_multipole_accept(ci->multipole, cj->multipole,
theta_crit_inv))
theta_crit_inv, 1))
scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci,
cj, 1);
}
......
......@@ -2525,43 +2525,24 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
* @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.
* @param Are we using the current value of CoM or the ones from the last
* rebuild ?
*/
__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));
}
/**
* @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];
double theta_crit_inv, int rebuild) {
const double r_crit_a =
(rebuild ? ma->r_max_rebuild : ma->r_max) * theta_crit_inv;
const double r_crit_b =
(rebuild ? mb->r_max_rebuild : mb->r_max) * theta_crit_inv;
const double dx = rebuild ? ma->CoM_rebuild[0] - mb->CoM_rebuild[0]
: ma->CoM[0] - mb->CoM[0];
const double dy = rebuild ? ma->CoM_rebuild[1] - mb->CoM_rebuild[1]
: ma->CoM[1] - mb->CoM[1];
const double dz = rebuild ? ma->CoM_rebuild[2] - mb->CoM_rebuild[2]
: ma->CoM[2] - mb->CoM[2];
const double r2 = dx * dx + dy * dy + dz * dz;
......
......@@ -458,7 +458,8 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
TIMER_TIC;
/* Can we use M-M interactions ? */
if (gravity_multipole_accept(ci->multipole, cj->multipole, theta_crit_inv)) {
if (gravity_multipole_accept(ci->multipole, cj->multipole, theta_crit_inv,
0)) {
/* MATTHIEU: make a symmetric M-M interaction function ! */
runner_dopair_grav_mm(r, ci, cj);
runner_dopair_grav_mm(r, cj, ci);
......@@ -638,11 +639,19 @@ 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_rebuild(ci->multipole, cj->multipole,
theta_crit_inv)) {
if (gravity_multipole_accept(ci->multipole, cj->multipole, theta_crit_inv,
0)) {
/* Go for a (non-symmetric) M-M calculation */
runner_dopair_grav_mm(r, ci, cj);
}
/* Is the criterion violated now but was OK at the last rebuild ? */
else if (gravity_multipole_accept(ci->multipole, cj->multipole,
theta_crit_inv, 1)) {
/* Alright, we have to take charge of that pair in a different way. */
runner_dopair_grav_mm(r, ci, cj);
}
}
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