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

Post-merge corrections

parent adb3d8d4
No related branches found
No related tags found
2 merge requests!212Gravity infrastructure,!172[WIP] Self gravity (Barnes-Hut version)
...@@ -1183,9 +1183,9 @@ void *runner_main(void *data) { ...@@ -1183,9 +1183,9 @@ void *runner_main(void *data) {
case task_type_grav_up: case task_type_grav_up:
runner_dograv_up(r, t->ci); runner_dograv_up(r, t->ci);
break; break;
case task_type_grav_down: /* case task_type_grav_down: */
runner_dograv_down(r, t->ci); /* runner_dograv_down(r, t->ci); */
break; /* break; */
case task_type_grav_external: case task_type_grav_external:
runner_do_grav_external(r, t->ci, 1); runner_do_grav_external(r, t->ci, 1);
break; break;
......
...@@ -25,178 +25,6 @@ ...@@ -25,178 +25,6 @@
#include "clocks.h" #include "clocks.h"
#include "part.h" #include "part.h"
/**
* @brief Compute the sorted gravity interactions between a cell pair.
*
* @param r The #runner.
* @param ci The first #cell.
* @param cj The second #cell.
*/
void runner_dopair_grav_new(struct runner *r, struct cell *ci,
struct cell *cj) {
struct engine *restrict e = r->e;
int pid, pjd, k, sid;
double rshift, shift[3] = {0.0, 0.0, 0.0}, nshift[3];
struct entry *restrict sort_i, *restrict sort_j;
struct gpart *restrict pi, *restrict pj, *restrict parts_i, *restrict parts_j;
double pix[3];
float dx[3], r2, h_max, di, dj;
int count_i, count_j, cnj, cnj_new;
const int ti_current = e->ti_current;
struct multipole m;
#ifdef VECTORIZE
int icount = 0;
float r2q[VEC_SIZE] __attribute__((aligned(16)));
float dxq[3 * VEC_SIZE] __attribute__((aligned(16)));
struct gpart *piq[VEC_SIZE], *pjq[VEC_SIZE];
#endif
TIMER_TIC
/* Anything to do here? */
if (ci->ti_end_min > ti_current && cj->ti_end_min > ti_current) return;
/* Get the sort ID. */
sid = space_getsid(e->s, &ci, &cj, shift);
/* Make sure the cells are sorted. */
runner_dogsort(r, ci, (1 << sid), 0);
runner_dogsort(r, cj, (1 << sid), 0);
/* Have the cells been sorted? */
if (!(ci->gsorted & (1 << sid)) || !(cj->gsorted & (1 << sid)))
error("Trying to interact unsorted cells.");
/* Get the cutoff shift. */
for (rshift = 0.0, k = 0; k < 3; k++)
rshift += shift[k] * runner_shift[3 * sid + k];
/* Pick-out the sorted lists. */
sort_i = &ci->gsort[sid * (ci->count + 1)];
sort_j = &cj->gsort[sid * (cj->count + 1)];
/* Get some other useful values. */
h_max =
sqrtf(ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2]) *
const_theta_max;
count_i = ci->gcount;
count_j = cj->gcount;
parts_i = ci->gparts;
parts_j = cj->gparts;
cnj = count_j;
multipole_reset(&m);
nshift[0] = -shift[0];
nshift[1] = -shift[1];
nshift[2] = -shift[2];
/* Loop over the parts in ci. */
for (pid = count_i - 1; pid >= 0; pid--) {
/* Get a hold of the ith part in ci. */
pi = &parts_i[sort_i[pid].i];
if (pi->ti_end > ti_current) continue;
di = sort_i[pid].d + h_max - rshift;
for (k = 0; k < 3; k++) pix[k] = pi->x[k] - shift[k];
/* Loop over the parts in cj. */
for (pjd = 0; pjd < cnj && sort_j[pjd].d < di; pjd++) {
/* Get a pointer to the jth particle. */
pj = &parts_j[sort_j[pjd].i];
/* Compute the pairwise distance. */
r2 = 0.0f;
for (k = 0; k < 3; k++) {
dx[k] = pix[k] - pj->x[k];
r2 += dx[k] * dx[k];
}
#ifndef VECTORIZE
// if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525 )
// message( "interacting particles pi=%lli and pj=%lli with r=%.3e in
// cells %lli/%lli." , pi->part->id , pj->part->id , sqrtf(r2) , ((long
// long int)ci) / sizeof(struct cell) , ((long long int)cj) /
// sizeof(struct cell) );
runner_iact_grav(r2, dx, pi, pj);
#else
/* Add this interaction to the queue. */
r2q[icount] = r2;
dxq[3 * icount + 0] = dx[0];
dxq[3 * icount + 1] = dx[1];
dxq[3 * icount + 2] = dx[2];
piq[icount] = pi;
pjq[icount] = pj;
icount += 1;
/* Flush? */
if (icount == VEC_SIZE) {
runner_iact_vec_grav(r2q, dxq, piq, pjq);
icount = 0;
}
#endif
} /* loop over the parts in cj. */
/* Set the new limit. */
cnj_new = pjd;
/* Add trailing parts to the multipole. */
for (pjd = cnj_new; pjd < cnj; pjd++) {
/* Add the part to the multipole. */
multipole_addpart(&m, &parts_j[sort_j[pjd].i]);
} /* add trailing parts to the multipole. */
/* Set the new cnj. */
cnj = cnj_new;
/* Interact the ith particle with the multipole. */
multipole_iact_mp(&m, pi, nshift);
} /* loop over the parts in ci. */
#ifdef VECTORIZE
/* Pick up any leftovers. */
if (icount > 0)
for (k = 0; k < icount; k++)
runner_iact_grav(r2q[k], &dxq[3 * k], piq[k], pjq[k]);
#endif
/* Re-set the multipole. */
multipole_reset(&m);
/* Loop over the parts in cj and interact with the multipole in ci. */
for (pid = count_i - 1, pjd = 0; pjd < count_j; pjd++) {
/* Get the position of pj along the axis. */
dj = sort_j[pjd].d - h_max + rshift;
/* Add any left-over parts in cell_i to the multipole. */
while (pid >= 0 && sort_i[pid].d < dj) {
/* Add this particle to the multipole. */
multipole_addpart(&m, &parts_i[sort_i[pid].i]);
/* Decrease pid. */
pid -= 1;
}
/* Interact pj with the multipole. */
multipole_iact_mp(&m, &parts_j[sort_j[pjd].i], shift);
} /* loop over the parts in cj and interact with the multipole. */
TIMER_TOC(TIMER_DOPAIR);
}
/** /**
* @brief Compute the recursive upward sweep, i.e. construct the * @brief Compute the recursive upward sweep, i.e. construct the
* multipoles in a cell hierarchy. * multipoles in a cell hierarchy.
...@@ -389,7 +217,7 @@ __attribute__((always_inline)) INLINE static void runner_dograv_pair_pp( ...@@ -389,7 +217,7 @@ __attribute__((always_inline)) INLINE static void runner_dograv_pair_pp(
* @brief Computes the interaction of all the particles in a cell directly * @brief Computes the interaction of all the particles in a cell directly
* *
* @param r The #runner. * @param r The #runner.
* @param ci The first #cell. * @param c The #cell.
* *
* @todo Use a local cache for the particles. * @todo Use a local cache for the particles.
*/ */
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment