Commit 5ea62658 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Updated all the calls to M2L to use the new symmetric M2L version

parent 94f1bef1
......@@ -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);
}
......@@ -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
......
......@@ -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;
}
}
......
......@@ -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,
......
......@@ -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);
......
......@@ -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
......
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