diff --git a/src/cell.c b/src/cell.c index 0702d0aa79c8bf97923eee955ef1216b7c2472a8..4454fe69376a71bdd2f6498837d2be54ebe3c5dd 100644 --- a/src/cell.c +++ b/src/cell.c @@ -2219,6 +2219,7 @@ void cell_make_multipoles(struct cell *c, integertime_t ti_current, gravity_reset(c->grav.multipole); if (c->split) { + /* Start by recursing */ for (int k = 0; k < 8; ++k) { if (c->progeny[k] != NULL) @@ -2281,9 +2282,18 @@ void cell_make_multipoles(struct cell *c, integertime_t ti_current, /* Take minimum of both limits */ c->grav.multipole->r_max = min(r_max, sqrt(dx * dx + dy * dy + dz * dz)); + /* We know the first-order multipole (dipole) is 0. */ + c->grav.multipole->m_pole.M_100 = 0.f; + c->grav.multipole->m_pole.M_010 = 0.f; + c->grav.multipole->m_pole.M_001 = 0.f; + + /* Compute the multipole power */ + gravity_multipole_compute_power(&c->grav.multipole->m_pole); + } else { if (c->grav.count > 0) { gravity_P2M(c->grav.multipole, c->grav.parts, c->grav.count, grav_props); + gravity_multipole_compute_power(&c->grav.multipole->m_pole); const double dx = c->grav.multipole->CoM[0] > c->loc[0] + c->width[0] * 0.5 ? c->grav.multipole->CoM[0] - c->loc[0] @@ -2377,6 +2387,7 @@ void cell_check_multipole(struct cell *c, if (c->grav.count > 0) { /* Brute-force calculation */ gravity_P2M(&ma, c->grav.parts, c->grav.count, grav_props); + gravity_multipole_compute_power(&ma.m_pole); /* Now compare the multipole expansion */ if (!gravity_multipole_equal(&ma, c->grav.multipole, tolerance)) { diff --git a/src/runner_doiact_grav.c b/src/runner_doiact_grav.c index 91abe9e91060b19fa6df1bbbb29e0de6e6442cad..81c091a8946f22a1557d29f09cb9ce1fb8c23d15 100644 --- a/src/runner_doiact_grav.c +++ b/src/runner_doiact_grav.c @@ -742,7 +742,7 @@ static INLINE void runner_dopair_grav_pm_truncated( * @param ci The first #cell. * @param cj The other #cell. * @param symmetric Are we updating both cells (1) or just ci (0) ? - * @param allow_mpole Are we allowing the use of P2M interactions ? + * @param allow_mpole Are we allowing the use of M2P interactions ? */ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj, const int symmetric, const int allow_mpole) { diff --git a/src/space.c b/src/space.c index 765d1d69decb1f59411642cd9a52cb29af32c678..51108737899c3ea476421af089843dfcb96b4e39 100644 --- a/src/space.c +++ b/src/space.c @@ -3664,6 +3664,9 @@ void space_split_recursive(struct space *s, struct cell *c, c->grav.multipole->m_pole.M_010 = 0.f; c->grav.multipole->m_pole.M_001 = 0.f; + /* Compute the multipole power */ + gravity_multipole_compute_power(&c->grav.multipole->m_pole); + } /* Deal with gravity */ } /* Split or let it be? */ @@ -3801,6 +3804,9 @@ void space_split_recursive(struct space *s, struct cell *c, gravity_P2M(c->grav.multipole, c->grav.parts, c->grav.count, e->gravity_properties); + /* Compute the multipole power */ + gravity_multipole_compute_power(&c->grav.multipole->m_pole); + } else { /* No gparts in that leaf cell */