diff --git a/src/runner.c b/src/runner.c index c59c80ffb048eb787a14e8096e6ace76377f5981..5921691fba4c20adb136fee5513ae381761a6c71 100644 --- a/src/runner.c +++ b/src/runner.c @@ -1183,9 +1183,9 @@ void *runner_main(void *data) { case task_type_grav_up: runner_dograv_up(r, t->ci); break; - case task_type_grav_down: - runner_dograv_down(r, t->ci); - break; + /* case task_type_grav_down: */ + /* runner_dograv_down(r, t->ci); */ + /* break; */ case task_type_grav_external: runner_do_grav_external(r, t->ci, 1); break; diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h index a4d34e802461753997c23108bfb5cefbac73f50f..5ea19adf5df3bcb714e61709e5984dd5da02000a 100644 --- a/src/runner_doiact_grav.h +++ b/src/runner_doiact_grav.h @@ -25,178 +25,6 @@ #include "clocks.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 * multipoles in a cell hierarchy. @@ -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 * * @param r The #runner. - * @param ci The first #cell. + * @param c The #cell. * * @todo Use a local cache for the particles. */