diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h index 3a8cc90557792375af6d21f90fe566b0b0660fc2..62f4f4a5be68924f8c1b159e0e20f50297190316 100644 --- a/src/gravity/Default/gravity_iact.h +++ b/src/gravity/Default/gravity_iact.h @@ -128,6 +128,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated( /** * @brief Computes the forces at a point generated by a multipole. * + * This assumes M_100 == M_010 == M_001 == 0. * This uses the quadrupole and trace of the octupole terms only and defaults to * the monopole if the code is compiled with low-order gravity only. * @@ -209,6 +210,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full( * @brief Computes the forces at a point generated by a multipole, truncated for * long-range periodicity. * + * This assumes M_100 == M_010 == M_001 == 0. * This uses the quadrupole term and trace of the octupole terms only and * defaults to the monopole if the code is compiled with low-order gravity only. * diff --git a/src/multipole.h b/src/multipole.h index 3e0533a6a6a2aade55c34e4d9c9d3f211fda0078..231a4f5f1646a7baf41766a431991782256e55d0 100644 --- a/src/multipole.h +++ b/src/multipole.h @@ -2263,9 +2263,6 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb, gp->num_interacted += lb->num_interacted; #endif - // MATTHIEU - // return; - /* Local accumulator */ double a_grav[3] = {0., 0., 0.}; double pot = 0.; @@ -2431,7 +2428,6 @@ __attribute__((always_inline)) INLINE static int gravity_M2P_accept( /* Multipole acceptance criterion (Dehnen 2002, eq.10) */ return (r2 * theta_crit2 * 0.1 > r_max2); - // return 0; } #endif /* SWIFT_MULTIPOLE_H */ diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h index 404b0a52df189f6f97a924f5767d8bb8e3b039a7..c8cf98070895253852acfd760af8b2077d77dbf4 100644 --- a/src/runner_doiact_grav.h +++ b/src/runner_doiact_grav.h @@ -30,9 +30,6 @@ #include "space_getsid.h" #include "timers.h" -static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci, - struct cell *cj, int symmetric); - /** * @brief Recursively propagate the multipoles down the tree by applying the * L2L and L2P kernels. @@ -47,10 +44,6 @@ static INLINE void runner_do_grav_down(struct runner *r, struct cell *c, /* Some constants */ const struct engine *e = r->e; - /* Cell properties */ - struct gpart *gparts = c->gparts; - const int gcount = c->gcount; - TIMER_TIC; #ifdef SWIFT_DEBUG_CHECKS @@ -103,6 +96,13 @@ static INLINE void runner_do_grav_down(struct runner *r, struct cell *c, if (!cell_are_gpart_drifted(c, e)) error("Un-drifted gparts"); + /* Cell properties */ + struct gpart *gparts = c->gparts; + const int gcount = c->gcount; + const struct grav_tensor *pot = &c->multipole->pot; + const double CoM[3] = {c->multipole->CoM[0], c->multipole->CoM[1], + c->multipole->CoM[2]}; + /* Apply accelerations to the particles */ for (int i = 0; i < gcount; ++i) { @@ -120,7 +120,7 @@ static INLINE void runner_do_grav_down(struct runner *r, struct cell *c, error("c->field tensor not initialised"); #endif /* Apply the kernel */ - gravity_L2P(&c->multipole->pot, c->multipole->CoM, gp); + gravity_L2P(pot, CoM, gp); } } } @@ -128,59 +128,6 @@ static INLINE void runner_do_grav_down(struct runner *r, struct cell *c, if (timer) TIMER_TOC(timer_dograv_down); } -/** - * @brief Computes the interaction of the field tensor in a cell with the - * multipole of another cell. - * - * @param r The #runner. - * @param ci The #cell with field tensor to interact. - * @param cj The #cell with the multipole. - */ -static INLINE void runner_dopair_grav_mm(struct runner *r, - struct cell *restrict ci, - struct cell *restrict cj) { - - /* Some constants */ - const struct engine *e = r->e; - const struct gravity_props *props = e->gravity_properties; - const int periodic = e->mesh->periodic; - const double dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]}; - const float r_s_inv = e->mesh->r_s_inv; - - TIMER_TIC; - - /* Anything to do here? */ - if (!cell_is_active_gravity(ci, e) || ci->nodeID != engine_rank) return; - - /* Short-cut to the multipole */ - const struct multipole *multi_j = &cj->multipole->m_pole; - -#ifdef SWIFT_DEBUG_CHECKS - if (ci == cj) error("Interacting a cell with itself using M2L"); - - if (multi_j->num_gpart == 0) - error("Multipole does not seem to have been set."); - - if (ci->multipole->pot.ti_init != e->ti_current) - error("ci->grav tensor not initialised."); -#endif - - /* Do we need to drift the multipole ? */ - if (cj->ti_old_multipole != e->ti_current) - error( - "Undrifted multipole cj->ti_old_multipole=%lld cj->nodeID=%d " - "ci->nodeID=%d e->ti_current=%lld", - cj->ti_old_multipole, cj->nodeID, ci->nodeID, e->ti_current); - - /* Let's interact at this level */ - gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM, - cj->multipole->CoM, props, periodic, dim, r_s_inv); - - runner_dopair_grav_pp(r, ci, cj, 0); - - TIMER_TOC(timer_dopair_grav_mm); -} - /** * @brief Compute the non-truncated gravity interactions between all particles * of a cell and the particles of the other cell. @@ -1160,6 +1107,60 @@ static INLINE void runner_doself_grav_pp(struct runner *r, struct cell *c) { TIMER_TOC(timer_doself_grav_pp); } +/** + * @brief Computes the interaction of the field tensor in a cell with the + * multipole of another cell. + * + * @param r The #runner. + * @param ci The #cell with field tensor to interact. + * @param cj The #cell with the multipole. + */ +static INLINE void runner_dopair_grav_mm(struct runner *r, + struct cell *restrict ci, + struct cell *restrict cj) { + + /* Some constants */ + const struct engine *e = r->e; + const struct gravity_props *props = e->gravity_properties; + const int periodic = e->mesh->periodic; + const double dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]}; + const float r_s_inv = e->mesh->r_s_inv; + + TIMER_TIC; + + /* Anything to do here? */ + if (!cell_is_active_gravity(ci, e) || ci->nodeID != engine_rank) return; + + /* Short-cut to the multipole */ + const struct multipole *multi_j = &cj->multipole->m_pole; + +#ifdef SWIFT_DEBUG_CHECKS + if (ci == cj) error("Interacting a cell with itself using M2L"); + + if (multi_j->num_gpart == 0) + error("Multipole does not seem to have been set."); + + if (ci->multipole->pot.ti_init != e->ti_current) + error("ci->grav tensor not initialised."); +#endif + + /* Do we need to drift the multipole ? */ + if (cj->ti_old_multipole != e->ti_current) + error( + "Undrifted multipole cj->ti_old_multipole=%lld cj->nodeID=%d " + "ci->nodeID=%d e->ti_current=%lld", + cj->ti_old_multipole, cj->nodeID, ci->nodeID, e->ti_current); + + /* Let's interact at this level */ + if (0) + gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM, + cj->multipole->CoM, props, periodic, dim, r_s_inv); + + runner_dopair_grav_pp(r, ci, cj, 0); + + TIMER_TOC(timer_dopair_grav_mm); +} + /** * @brief Computes the interaction of all the particles in a cell with all the * particles of another cell.