diff --git a/src/cell.c b/src/cell.c
index 4d9f25a4e286355c58710271ad0f938857873154..fe2f78d078993448474ec0c01db902dc512908bd 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -6382,6 +6382,8 @@ void cell_reorder_extra_gparts(struct cell *c, struct part *parts,
/**
* @brief Can we use the MM interactions fo a given pair of cells?
*
+ * This uses the information from the last tree-build.
+ *
* @param ci The first #cell.
* @param cj The second #cell.
* @param e The #engine.
@@ -6390,8 +6392,7 @@ void cell_reorder_extra_gparts(struct cell *c, struct part *parts,
int cell_can_use_pair_mm(const struct cell *ci, const struct cell *cj,
const struct engine *e, const struct space *s) {
- const double theta_crit = e->gravity_properties->theta_crit;
- const double theta_crit2 = theta_crit * theta_crit;
+ const struct gravity_props *props = e->gravity_properties;
const int periodic = s->periodic;
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
@@ -6412,76 +6413,6 @@ int cell_can_use_pair_mm(const struct cell *ci, const struct cell *cj,
}
const double r2 = dx * dx + dy * dy + dz * dz;
- const double epsilon_i = multi_i->m_pole.max_softening;
- const double epsilon_j = multi_j->m_pole.max_softening;
-
- return gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2,
- epsilon_i, epsilon_j);
-}
-
-/**
- * @brief Can we use the MM interactions fo a given pair of cells?
- *
- * This function uses the information gathered in the multipole at rebuild
- * time and not the current position and radius of the multipole.
- *
- * @param ci The first #cell.
- * @param cj The second #cell.
- * @param e The #engine.
- * @param s The #space.
- */
-int cell_can_use_pair_mm_rebuild(const struct cell *ci, const struct cell *cj,
- const struct engine *e,
- const struct space *s) {
- const double theta_crit = e->gravity_properties->theta_crit;
- const double theta_crit2 = theta_crit * theta_crit;
- const int periodic = s->periodic;
- const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
-
- /* Recover the multipole information */
- const struct gravity_tensors *const multi_i = ci->grav.multipole;
- const struct gravity_tensors *const multi_j = cj->grav.multipole;
-
-#ifdef SWIFT_DEBUG_CHECKS
-
- if (multi_i->CoM_rebuild[0] < ci->loc[0] ||
- multi_i->CoM_rebuild[0] > ci->loc[0] + ci->width[0])
- error("Invalid multipole position ci");
- if (multi_i->CoM_rebuild[1] < ci->loc[1] ||
- multi_i->CoM_rebuild[1] > ci->loc[1] + ci->width[1])
- error("Invalid multipole position ci");
- if (multi_i->CoM_rebuild[2] < ci->loc[2] ||
- multi_i->CoM_rebuild[2] > ci->loc[2] + ci->width[2])
- error("Invalid multipole position ci");
-
- if (multi_j->CoM_rebuild[0] < cj->loc[0] ||
- multi_j->CoM_rebuild[0] > cj->loc[0] + cj->width[0])
- error("Invalid multipole position cj");
- if (multi_j->CoM_rebuild[1] < cj->loc[1] ||
- multi_j->CoM_rebuild[1] > cj->loc[1] + cj->width[1])
- error("Invalid multipole position cj");
- if (multi_j->CoM_rebuild[2] < cj->loc[2] ||
- multi_j->CoM_rebuild[2] > cj->loc[2] + cj->width[2])
- error("Invalid multipole position cj");
-
-#endif
-
- /* Get the distance between the CoMs */
- double dx = multi_i->CoM_rebuild[0] - multi_j->CoM_rebuild[0];
- double dy = multi_i->CoM_rebuild[1] - multi_j->CoM_rebuild[1];
- double dz = multi_i->CoM_rebuild[2] - multi_j->CoM_rebuild[2];
-
- /* Apply BC */
- if (periodic) {
- dx = nearest(dx, dim[0]);
- dy = nearest(dy, dim[1]);
- dz = nearest(dz, dim[2]);
- }
- const double r2 = dx * dx + dy * dy + dz * dz;
-
- const double epsilon_i = multi_i->m_pole.max_softening;
- const double epsilon_j = multi_j->m_pole.max_softening;
-
- return gravity_M2L_accept(multi_i->r_max_rebuild, multi_j->r_max_rebuild,
- theta_crit2, r2, epsilon_i, epsilon_j);
+ return gravity_M2L_accept_symmetric(props, multi_i, multi_j, r2,
+ /*use_rebuild_sizes=*/1);
}
diff --git a/src/cell.h b/src/cell.h
index 79e4d875f3da38a08143551ca90e494a820d24b6..a78e69fea5b3e2f76a86ed5a52b76644285974e7 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -960,8 +960,6 @@ void cell_reorder_extra_gparts(struct cell *c, struct part *parts,
void cell_reorder_extra_sparts(struct cell *c, const ptrdiff_t sparts_offset);
int cell_can_use_pair_mm(const struct cell *ci, const struct cell *cj,
const struct engine *e, const struct space *s);
-int cell_can_use_pair_mm_rebuild(const struct cell *ci, const struct cell *cj,
- const struct engine *e, const struct space *s);
/**
* @brief Compute the square of the minimal distance between any two points in
diff --git a/src/engine.c b/src/engine.c
index d1760d0b3fe52d30c0f8b8c9aa38eced5ec09d89..1dc486e2da220885b3d1e514d0ae2515333aa8b8 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3100,8 +3100,6 @@ void engine_makeproxies(struct engine *e) {
const int with_hydro = (e->policy & engine_policy_hydro);
const int with_gravity = (e->policy & engine_policy_self_gravity);
const double theta_crit_inv = 1. / e->gravity_properties->theta_crit;
- const double theta_crit2 =
- e->gravity_properties->theta_crit * e->gravity_properties->theta_crit;
const double max_mesh_dist = e->mesh->r_cut_max;
const double max_mesh_dist2 = max_mesh_dist * max_mesh_dist;
@@ -3245,26 +3243,29 @@ void engine_makeproxies(struct engine *e) {
sqrt(min_dist_centres2) - 2. * delta_CoM;
const double min_dist_CoM2 = min_dist_CoM * min_dist_CoM;
- /* We also assume that the softening is negligible compared
- to the cell size */
- const double epsilon_i = 0.;
- const double epsilon_j = 0.;
+ error("MATTHIEU: Need to implement a good MAC for proxies!");
+
+ /* /\* We also assume that the softening is negligible
+ * compared */
+ /* to the cell size *\/ */
+ /* const double epsilon_i = 0.; */
+ /* const double epsilon_j = 0.; */
/* Are we beyond the distance where the truncated forces are 0
* but not too far such that M2L can be used? */
if (periodic) {
- if ((min_dist_CoM2 < max_mesh_dist2) &&
- (!gravity_M2L_accept(r_max, r_max, theta_crit2,
- min_dist_CoM2, epsilon_i,
- epsilon_j)))
+ if ((min_dist_CoM2 < max_mesh_dist2)) // &&
+ /* (!gravity_M2L_accept(r_max, r_max, theta_crit2, */
+ /* min_dist_CoM2, epsilon_i, */
+ /* epsilon_j))) */
proxy_type |= (int)proxy_cell_type_gravity;
} else {
- if (!gravity_M2L_accept(r_max, r_max, theta_crit2,
- min_dist_CoM2, epsilon_i,
- epsilon_j))
+ if (1 /* !gravity_M2L_accept(r_max, r_max, theta_crit2, */
+ /* min_dist_CoM2, epsilon_i, */
+ /* epsilon_j) */)
proxy_type |= (int)proxy_cell_type_gravity;
}
}
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index 9e48e5826fc13256647235ffa1c5dad00adb76bb..73f1d6d44d4d5a1afc037af16502afb61f25ab2e 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -1381,7 +1381,7 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
if (periodic && min_radius2 > max_distance2) continue;
/* Are the cells too close for a MM interaction ? */
- if (!cell_can_use_pair_mm_rebuild(ci, cj, e, s)) {
+ if (!cell_can_use_pair_mm(ci, cj, e, s)) {
/* Ok, we need to add a direct pair calculation */
scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0,
diff --git a/src/runner_doiact_grav.c b/src/runner_doiact_grav.c
index b36bedd09042f0c6df07a89e92e30113c2cdb014..5a2f18f18f0c2bf2909ddc26a227d2d287994b56 100644
--- a/src/runner_doiact_grav.c
+++ b/src/runner_doiact_grav.c
@@ -538,19 +538,7 @@ static INLINE void runner_dopair_grav_pm_full(
const float r2 = dx * dx + dy * dy + dz * dz;
-#ifdef SWIFT_DEBUG_CHECKS
- const float r_max_j = cj->grav.multipole->r_max;
- const float r_max2 = r_max_j * r_max_j;
- const double theta_crit = e->gravity_properties->theta_crit;
- const double theta_crit2 = theta_crit * theta_crit;
-
- /* Note: 0.99 and 1.1 to avoid FP rounding false-positives */
- if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.1, r2, 0.99 * h_i))
- error(
- "use_mpole[i] set when M2P accept fails CoM=[%e %e %e] pos=[%e %e "
- "%e], rmax=%e r=%e epsilon=%e",
- CoM_j[0], CoM_j[1], CoM_j[2], x_i, y_i, z_i, r_max_j, sqrtf(r2), h_i);
-#endif
+ // MATTHIEU: Do we want a check that the particle can indeed use M2P?
/* Interact! */
float f_x, f_y, f_z, pot_ij;
@@ -682,19 +670,7 @@ static INLINE void runner_dopair_grav_pm_truncated(
const float r2 = dx * dx + dy * dy + dz * dz;
-#ifdef SWIFT_DEBUG_CHECKS
- const float r_max_j = cj->grav.multipole->r_max;
- const float r_max2 = r_max_j * r_max_j;
- const double theta_crit = e->gravity_properties->theta_crit;
- const double theta_crit2 = theta_crit * theta_crit;
-
- /* 0.99 and 1.1 to avoid FP rounding false-positives */
- if (!gravity_M2P_accept(r_max2, theta_crit2 * 1.1, r2, 0.99 * h_i))
- error(
- "use_mpole[i] set when M2P accept fails CoM=[%e %e %e] pos=[%e %e "
- "%e], rmax=%e",
- CoM_j[0], CoM_j[1], CoM_j[2], x_i, y_i, z_i, r_max_j);
-#endif
+ // MATTHIEU: Do we want a check that the particle can indeed use M2P?
/* Interact! */
float f_x, f_y, f_z, pot_ij;
@@ -791,8 +767,6 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj,
/* Recover the multipole info and shift the CoM locations */
const float rmax_i = ci->grav.multipole->r_max;
const float rmax_j = cj->grav.multipole->r_max;
- const float rmax2_i = rmax_i * rmax_i;
- const float rmax2_j = rmax_j * rmax_j;
const struct multipole *multi_i = &ci->grav.multipole->m_pole;
const struct multipole *multi_j = &cj->grav.multipole->m_pole;
const float CoM_i[3] = {(float)(ci->grav.multipole->CoM[0] - shift_i[0]),
@@ -820,10 +794,12 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj,
/* Fill the caches */
gravity_cache_populate(e->max_active_bin, allow_mpole, periodic, dim,
ci_cache, ci->grav.parts, gcount_i, gcount_padded_i,
- shift_i, CoM_j, rmax2_j, ci, e->gravity_properties);
+ shift_i, CoM_j, cj->grav.multipole, ci,
+ e->gravity_properties);
gravity_cache_populate(e->max_active_bin, allow_mpole, periodic, dim,
cj_cache, cj->grav.parts, gcount_j, gcount_padded_j,
- shift_j, CoM_i, rmax2_i, cj, e->gravity_properties);
+ shift_j, CoM_i, ci->grav.multipole, cj,
+ e->gravity_properties);
/* Can we use the Newtonian version or do we need the truncated one ? */
if (!periodic) {
@@ -1573,7 +1549,6 @@ void runner_dopair_recursive_grav_pm(struct runner *r, struct cell *ci,
/* Recover the multipole info and the CoM locations */
const struct multipole *multi_j = &cj->grav.multipole->m_pole;
- const float r_max = cj->grav.multipole->r_max;
const float CoM_j[3] = {(float)(cj->grav.multipole->CoM[0]),
(float)(cj->grav.multipole->CoM[1]),
(float)(cj->grav.multipole->CoM[2])};
@@ -1581,7 +1556,7 @@ void runner_dopair_recursive_grav_pm(struct runner *r, struct cell *ci,
/* Fill the cache */
gravity_cache_populate_all_mpole(
e->max_active_bin, periodic, dim, ci_cache, ci->grav.parts, gcount_i,
- gcount_padded_i, ci, CoM_j, r_max * r_max, e->gravity_properties);
+ gcount_padded_i, ci, CoM_j, cj->grav.multipole, e->gravity_properties);
/* Can we use the Newtonian version or do we need the truncated one ? */
if (!periodic) {
@@ -1618,15 +1593,13 @@ void runner_dopair_recursive_grav_pm(struct runner *r, struct cell *ci,
* @param gettimer Are we timing this ?
*/
void runner_dopair_recursive_grav(struct runner *r, struct cell *ci,
- struct cell *cj, int gettimer) {
+ struct cell *cj, const int gettimer) {
/* Some constants */
const struct engine *e = r->e;
const int nodeID = e->nodeID;
const int periodic = e->mesh->periodic;
const double dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]};
- const double theta_crit = e->gravity_properties->theta_crit;
- const double theta_crit2 = theta_crit * theta_crit;
const double max_distance = e->mesh->r_cut_max;
/* Anything to do here? */
@@ -1704,9 +1677,8 @@ void runner_dopair_recursive_grav(struct runner *r, struct cell *ci,
* option... */
/* Can we use M-M interactions ? */
- if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2,
- multi_i->m_pole.max_softening,
- multi_j->m_pole.max_softening)) {
+ if (gravity_M2L_accept_symmetric(e->gravity_properties, multi_i, multi_j, r2,
+ /* use_rebuild_sizes=*/0)) {
/* Go M-M */
runner_dopair_grav_mm(r, ci, cj);
@@ -1786,7 +1758,7 @@ void runner_dopair_recursive_grav(struct runner *r, struct cell *ci,
* @param gettimer Are we timing this ?
*/
void runner_doself_recursive_grav(struct runner *r, struct cell *c,
- int gettimer) {
+ const int gettimer) {
/* Some constants */
const struct engine *e = r->e;
@@ -1837,14 +1809,13 @@ void runner_doself_recursive_grav(struct runner *r, struct cell *c,
* @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) {
+void runner_do_grav_long_range(struct runner *r, struct cell *ci,
+ const int timer) {
/* Some constants */
const struct engine *e = r->e;
const int periodic = e->mesh->periodic;
const double dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]};
- const double theta_crit = e->gravity_properties->theta_crit;
- const double theta_crit2 = theta_crit * theta_crit;
const double max_distance2 = e->mesh->r_cut_max * e->mesh->r_cut_max;
TIMER_TIC;
@@ -1934,10 +1905,8 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
const double r2_rebuild = dx_r * dx_r + dy_r * dy_r + dz_r * dz_r;
/* Are we in charge of this cell pair? */
- if (gravity_M2L_accept(multi_top->r_max_rebuild, multi_j->r_max_rebuild,
- theta_crit2, r2_rebuild,
- multi_top->m_pole.max_softening,
- multi_j->m_pole.max_softening)) {
+ if (gravity_M2L_accept_symmetric(e->gravity_properties, multi_top, multi_j,
+ r2_rebuild, /*use_rebuild_sizes=*/1)) {
/* Call the PM interaction fucntion on the active sub-cells of ci */
runner_dopair_grav_mm_nonsym(r, ci, cj);
diff --git a/src/scheduler.c b/src/scheduler.c
index df6f908a3f620d9d6abb3c9e07d179d2d30b1f96..268ee5023d11e0acb317a084987ccfaff4647fff 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -872,8 +872,9 @@ static void scheduler_splittask_gravity(struct task *t, struct scheduler *s) {
for (int j = 0; j < 8; j++) {
if (cj->progeny[j] != NULL) {
/* Can we use a M-M interaction here? */
- if (cell_can_use_pair_mm_rebuild(ci->progeny[i],
- cj->progeny[j], e, sp)) {
+ if (cell_can_use_pair_mm(ci->progeny[i], cj->progeny[j], e,
+ sp)) {
+
/* Flag this pair as being treated by the M-M task.
* We use the 64 bits in the task->flags field to store
* this information. The corresponding taks will unpack