Skip to content
Snippets Groups Projects
Commit 5026f624 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Restored the long-range gravity task. Uses the normal MAC again.

parent cf17832d
No related branches found
No related tags found
1 merge request!433Correct wrapping of multipoles in FFT task
Hide whitespace changes
Inline Side-by-side
src/runner_doiact_grav.h
+ 40
− 84
View file @ 5026f624
......@@ -1145,24 +1145,12 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
const struct gravity_props *props = e->gravity_properties;
const int periodic = s->periodic;
const double cell_width = s->width[0];
const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
/* const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]}; */
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const double theta_crit2 = props->theta_crit2;
const double max_distance = props->a_smooth * props->r_cut_max * cell_width;
const double max_distance2 = max_distance * max_distance;
#if (ICHECK != 0)
int check = 0;
int direct_ngbs = 0;
int direct_ngbs_gpart = 0;
int other_ngbs_gpart = 0;
for (int i = 0; i < ci->gcount; ++i)
if (ci->gparts[i].id_or_neg_offset == ICHECK) {
message("Found gpart");
check = 1;
}
#endif
#ifdef SWIFT_DEBUG_CHECKS
long long counter = 0;
#endif
......@@ -1186,15 +1174,15 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
/* Recover the local multipole */
struct gravity_tensors *const multi_i = ci->multipole;
const double CoM_i[3] = {multi_i->CoM[0], multi_i->CoM[1], multi_i->CoM[2]};
/* const double CoM_rebuild_i[3] = {multi_i->CoM_rebuild[0], */
/* multi_i->CoM_rebuild[1], */
/* multi_i->CoM_rebuild[2]}; */
const double CoM_rebuild_i[3] = {multi_i->CoM_rebuild[0],
multi_i->CoM_rebuild[1],
multi_i->CoM_rebuild[2]};
/* Get the cell index. MATTHIEU */
const int cid = (ci - cells); // / sizeof(struct cell);
const int i = cid / (cdim[1] * cdim[2]);
const int j = (cid / cdim[2]) % cdim[1];
const int k = cid % cdim[2];
/* /\* Get the cell index. MATTHIEU *\/ */
/* const int cid = (ci - cells); // / sizeof(struct cell); */
/* const int i = cid / (cdim[1] * cdim[2]); */
/* const int j = (cid / cdim[2]) % cdim[1]; */
/* const int k = cid % cdim[2]; */
/* Loop over all the top-level cells and go for a M-M interaction if
* well-separated */
......@@ -1211,53 +1199,36 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
counter += multi_j->m_pole.num_gpart;
#endif
// MATTHIEU
const int cjd = (cj - cells); // / sizeof(struct cell);
const int ii = cjd / (cdim[1] * cdim[2]);
const int jj = (cjd / cdim[2]) % cdim[1];
const int kk = cjd % cdim[2];
/* if(check) message("cid=%d cjd=%d cj->nodeID=%d cj->gcount=%d", */
/* cid, cjd, cj->nodeID, cj->gcount); */
/* /\* Get the distance between the CoMs at the last rebuild*\/ */
/* double dx_r = CoM_rebuild_i[0] - multi_j->CoM_rebuild[0]; */
/* double dy_r = CoM_rebuild_i[1] - multi_j->CoM_rebuild[1]; */
/* double dz_r = CoM_rebuild_i[2] - multi_j->CoM_rebuild[2]; */
/* /\* Apply BC *\/ */
/* if (periodic) { */
/* dx_r = nearest(dx_r, dim[0]); */
/* dy_r = nearest(dy_r, dim[1]); */
/* dz_r = nearest(dz_r, dim[2]); */
/* } */
/* const double r2_rebuild = dx_r * dx_r + dy_r * dy_r + dz_r * dz_r; */
/* Are we in charge of this cell pair? MATTHIEU*/
/* if (gravity_M2L_accept(multi_i->r_max_rebuild, multi_j->r_max_rebuild, */
/* theta_crit2, r2_rebuild)) { */
if ((abs(i - ii) <= 1 || abs(i - ii - cdim[0]) <= 1 ||
abs(i - ii + cdim[0]) <= 1) &&
(abs(j - jj) <= 1 || abs(j - jj - cdim[1]) <= 1 ||
abs(j - jj + cdim[1]) <= 1) &&
(abs(k - kk) <= 1 || abs(k - kk - cdim[2]) <= 1 ||
abs(k - kk + cdim[2]) <= 1)) {
#if (ICHECK != 0)
if (check) {
++direct_ngbs;
direct_ngbs_gpart += cj->multipole->m_pole.num_gpart;
message(
"Found direct neighbour %d: (i,j,k)=(%d,%d,%d) "
"(ii,jj,kk)=(%d,%d,%d) nodeID=%d",
direct_ngbs, i, j, k, ii, jj, kk, cj->nodeID);
}
#endif
} else {
#if (ICHECK != 0)
if (check) other_ngbs_gpart += cj->multipole->m_pole.num_gpart;
#endif
/* // MATTHIEU */
/* const int cjd = (cj - cells); // / sizeof(struct cell); */
/* const int ii = cjd / (cdim[1] * cdim[2]); */
/* const int jj = (cjd / cdim[2]) % cdim[1]; */
/* const int kk = cjd % cdim[2]; */
/* /\* if(check) message("cid=%d cjd=%d cj->nodeID=%d cj->gcount=%d", *\/ */
/* /\* cid, cjd, cj->nodeID, cj->gcount); *\/ */
/* Get the distance between the CoMs at the last rebuild*/
double dx_r = CoM_rebuild_i[0] - multi_j->CoM_rebuild[0];
double dy_r = CoM_rebuild_i[1] - multi_j->CoM_rebuild[1];
double dz_r = CoM_rebuild_i[2] - multi_j->CoM_rebuild[2];
/* Apply BC */
if (periodic) {
dx_r = nearest(dx_r, dim[0]);
dy_r = nearest(dy_r, dim[1]);
dz_r = nearest(dz_r, dim[2]);
}
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_i->r_max_rebuild, multi_j->r_max_rebuild,
theta_crit2, r2_rebuild)) {
/* if ((abs(i - ii) <= 1 || abs(i - ii - cdim[0]) <= 1 || */
/* abs(i - ii + cdim[0]) <= 1) && */
/* (abs(j - jj) <= 1 || abs(j - jj - cdim[1]) <= 1 || */
/* abs(j - jj + cdim[1]) <= 1) && */
/* (abs(k - kk) <= 1 || abs(k - kk - cdim[2]) <= 1 || */
/* abs(k - kk + cdim[2]) <= 1)) { */
/* Let's compute the current distance between the cell pair*/
double dx = CoM_i[0] - multi_j->CoM[0];
......@@ -1281,7 +1252,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
#endif
continue;
}
/* Check the multipole acceptance criterion */
if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2)) {
......@@ -1295,21 +1266,6 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
} /* We are in charge of this pair */
} /* Loop over top-level cells */
#ifdef SWIFT_DEBUG_CHECKS
counter += ci->multipole->m_pole.num_gpart;
if (counter != e->total_nr_gparts)
error("Not found the right number of particles in top-level interactions");
#endif
#ifdef ICHECK
if (check)
message(
"Interacted with %d indirectly and ignored %d direct interactions "
"(counter=%lld) nr_cells=%d total=%lld",
other_ngbs_gpart, direct_ngbs_gpart, counter, nr_cells,
e->total_nr_gparts);
#endif
if (timer) TIMER_TOC(timer_dograv_long_range);
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment