diff --git a/src/engine.c b/src/engine.c
index 04124ec5aa8435b2551a19de999f37886020385e..2541b3dbde1b84b58fcc2bab54aa6d0b632ab820 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2990,27 +2990,16 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
}
#ifdef SWIFT_DEBUG_CHECKS
- /* Let's store the total mass in the system for future checks */
- e->s->total_mass = 0.;
- for (size_t i = 0; i < s->nr_gparts; ++i)
- e->s->total_mass += s->gparts[i].mass;
-#ifdef WITH_MPI
- if (MPI_Allreduce(MPI_IN_PLACE, &e->s->total_mass, 1, MPI_DOUBLE, MPI_SUM,
- MPI_COMM_WORLD) != MPI_SUCCESS)
- error("Failed to all-reduce total mass in the system.");
-#endif
- if (e->nodeID == 0) message("Total mass in the system: %e", e->s->total_mass);
-
/* Check that we have the correct total mass in the top-level multipoles */
+ size_t num_gpart_mpole = 0;
if (e->policy & engine_policy_self_gravity) {
- double mass = 0.;
for (int i = 0; i < e->s->nr_cells; ++i)
- mass += e->s->cells_top[i].multipole->m_pole.M_000;
- if (fabs(mass - e->s->total_mass) > e->s->total_mass / e->s->nr_gparts)
+ num_gpart_mpole += e->s->cells_top[i].multipole->m_pole.num_gpart;
+ if (num_gpart_mpole != e->s->nr_gparts)
error(
- "Total mass in multipoles does not match particle content. part=%e "
- "m-poles=%e",
- e->s->total_mass, mass);
+ "Multipoles don't contain the total number of gpart s->nr_gpart=%zd, "
+ "m_poles=%zd",
+ e->s->nr_gparts, num_gpart_mpole);
}
#endif
@@ -3101,15 +3090,12 @@ void engine_step(struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check that we have the correct total mass in the top-level multipoles */
+ size_t num_gpart_mpole = 0;
if (e->policy & engine_policy_self_gravity) {
- double mass = 0.;
for (int i = 0; i < e->s->nr_cells; ++i)
- mass += e->s->cells_top[i].multipole->m_pole.M_000;
- if (fabs(mass - e->s->total_mass) > e->s->total_mass / e->s->nr_gparts)
- error(
- "Total mass in multipoles does not match particle content. part=%e "
- "m-poles=%e",
- e->s->total_mass, mass);
+ num_gpart_mpole += e->s->cells_top[i].multipole->m_pole.num_gpart;
+ if (num_gpart_mpole != e->s->nr_gparts)
+ error("Multipoles don't contain the total number of gpart");
}
#endif
diff --git a/src/gravity.c b/src/gravity.c
index ee7b9afa02f106e0036f97388ad0bc5695da217e..e91773038c91527615b82051e7e90b85e51903cb 100644
--- a/src/gravity.c
+++ b/src/gravity.c
@@ -149,10 +149,10 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
fabsf(err_rel[2]) > rel_tol)
error(
"Error too large ! gp->a_grav=[%e %e %e] gp->a_exact=[%e %e %e], "
- "gp->mass_interacted=%e",
+ "gp->num_interacted=%lld",
gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2],
gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2],
- gpi->mass_interacted);
+ gpi->num_interacted);
/* Construct some statistics */
for (int k = 0; k < 3; ++k) {
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 312a9b61b1caa97b9ccd9e2c092ca406ab374da2..88ccd7159284d4b479fd8ff14bcc01bbc0383f06 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -62,7 +62,7 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
gp->a_grav[2] = 0.f;
#ifdef SWIFT_DEBUG_CHECKS
- gp->mass_interacted = 0.;
+ gp->num_interacted = 0;
#endif
}
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 00ae0f5b05cd95750c34b60e2353a9fc1d0a5c32..785ea4e86866b645b8ff72491adfc8c821f683af 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -52,8 +52,8 @@ struct gpart {
#ifdef SWIFT_DEBUG_CHECKS
- /* Total mass this gpart interacted with */
- double mass_interacted;
+ /* Numer of gparts this gpart interacted with */
+ long long num_interacted;
/* Time of the last drift */
integertime_t ti_drift;
diff --git a/src/multipole.h b/src/multipole.h
index 8dc3af9610432374bfa67d173ece7b9b1cbb931c..399e058c2efd7f69d8ea5941794b0df76a5388b4 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -69,8 +69,8 @@ struct grav_tensor {
#endif
#ifdef SWIFT_DEBUG_CHECKS
- /* Total mass this gpart interacted with */
- double mass_interacted;
+ /* Total number of gpart this field tensor interacted with */
+ long long num_interacted;
#endif
};
@@ -106,6 +106,13 @@ struct multipole {
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#error "Missing implementation for order >3"
#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Total number of gpart in this multipole */
+ long long num_gpart;
+
+#endif
};
/**
@@ -169,8 +176,8 @@ INLINE static void gravity_field_tensors_init(struct grav_tensor *l) {
INLINE static void gravity_field_tensors_add(struct grav_tensor *la,
const struct grav_tensor *lb) {
#ifdef SWIFT_DEBUG_CHECKS
- if (lb->mass_interacted == 0.f) error("Adding tensors that did not interact");
- la->mass_interacted += lb->mass_interacted;
+ if (lb->num_interacted == 0) error("Adding tensors that did not interact");
+ la->num_interacted += lb->num_interacted;
#endif
/* Add 0th order terms */
@@ -336,6 +343,10 @@ INLINE static void gravity_multipole_add(struct multipole *ma,
ma->M_012 += mb->M_012;
ma->M_111 += mb->M_111;
#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+ ma->num_gpart += mb->num_gpart;
+#endif
}
/**
@@ -648,6 +659,10 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#error "Missing implementation for order >3"
#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+ m->m_pole.num_gpart = gcount;
+#endif
}
/**
@@ -740,6 +755,10 @@ INLINE static void gravity_M2M(struct multipole *m_a,
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#error "Missing implementation for order >3"
#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+ m_a->num_gpart = m_b->num_gpart;
+#endif
}
/**
@@ -888,7 +907,8 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
#endif
#ifdef SWIFT_DEBUG_CHECKS
- l_b->mass_interacted += m_a->M_000;
+
+ l_b->num_interacted += m_a->num_gpart;
#endif
}
@@ -1010,9 +1030,8 @@ INLINE static void gravity_L2L(struct grav_tensor *la,
#endif
#ifdef SWIFT_DEBUG_CHECKS
- if (lb->mass_interacted == 0.f)
- error("Shifting tensors that did not interact");
- la->mass_interacted = lb->mass_interacted;
+ if (lb->num_interacted == 0) error("Shifting tensors that did not interact");
+ la->num_interacted = lb->num_interacted;
#endif
}
@@ -1020,21 +1039,22 @@ INLINE static void gravity_L2L(struct grav_tensor *la,
* @brief Applies the #grav_tensor to a #gpart.
*
* Corresponds to equation (28a).
+ *
+ * @param lb The gravity field tensor to apply.
+ * @param loc The position of the gravity field tensor.
+ * @param gp The #gpart to update.
*/
-INLINE static void gravity_L2P(const struct gravity_tensors *l_b,
- struct gpart *gp) {
+INLINE static void gravity_L2P(const struct grav_tensor *lb,
+ const double loc[3], struct gpart *gp) {
#ifdef SWIFT_DEBUG_CHECKS
- if (l_b->pot.mass_interacted == 0.f)
- error("Interacting with empty field tensor");
- gp->mass_interacted += l_b->pot.mass_interacted;
+ if (lb->num_interacted == 0) error("Interacting with empty field tensor");
+ gp->num_interacted += lb->num_interacted;
#endif
- const struct grav_tensor *lb = &l_b->pot;
-
/* Distance to the multipole */
- const double dx[3] = {gp->x[0] - l_b->CoM[0], gp->x[1] - l_b->CoM[1],
- gp->x[2] - l_b->CoM[2]};
+ const double dx[3] = {gp->x[0] - loc[0], gp->x[1] - loc[1],
+ gp->x[2] - loc[2]};
/* 0th order interaction */
gp->a_grav[0] += X_000(dx) * lb->F_100;
diff --git a/src/runner.c b/src/runner.c
index f6480ff6173729d57980cd55d642d50d6fc2e1ad..4d8c839549e3057c3a4128626a635f33ee77650e 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1372,12 +1372,12 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
#ifdef SWIFT_DEBUG_CHECKS
if (e->policy & engine_policy_self_gravity) {
- gp->mass_interacted += gp->mass;
- if (fabs(gp->mass_interacted - e->s->total_mass) > gp->mass)
+ gp->num_interacted++;
+ if (gp->num_interacted - e->s->nr_gparts)
error(
"g-particle did not interact gravitationally with all other "
- "particles gp->mass_interacted=%e, total_mass=%e, gp->mass=%e",
- gp->mass_interacted, e->s->total_mass, gp->mass);
+ "particles gp->num_interacted=%lld, total_gparts=%zd",
+ gp->num_interacted, e->s->nr_gparts);
}
#endif
}
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index ebca99f7f838e14018bbdb87ea68b8efe73bc350..4ff2abcfb095965518c185446b2703816dccbbc2 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -60,7 +60,8 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
/* Apply accelerations to the particles */
for (int i = 0; i < gcount; ++i) {
struct gpart *gp = &gparts[i];
- if (gpart_is_active(gp, e)) gravity_L2P(c->multipole, gp);
+ if (gpart_is_active(gp, e))
+ gravity_L2P(&c->multipole->pot, c->multipole->CoM, gp);
}
}
@@ -250,7 +251,7 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
#ifdef SWIFT_DEBUG_CHECKS
- gpi->mass_interacted += gpj->mass;
+ gpi->num_interacted++;
#endif
}
}
@@ -279,7 +280,7 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
#ifdef SWIFT_DEBUG_CHECKS
- gpj->mass_interacted += gpi->mass;
+ gpj->num_interacted++;
#endif
}
}
@@ -350,8 +351,8 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
#ifdef SWIFT_DEBUG_CHECKS
- gpi->mass_interacted += gpj->mass;
- gpj->mass_interacted += gpi->mass;
+ gpi->num_interacted++;
+ gpj->num_interacted++;
#endif
} else {
@@ -361,7 +362,7 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
#ifdef SWIFT_DEBUG_CHECKS
- gpi->mass_interacted += gpj->mass;
+ gpi->num_interacted++;
#endif
} else if (gpart_is_active(gpj, e)) {
@@ -372,7 +373,7 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
#ifdef SWIFT_DEBUG_CHECKS
- gpj->mass_interacted += gpi->mass;
+ gpj->num_interacted++;
#endif
}
}
diff --git a/src/space.h b/src/space.h
index e3886364bfa7b2f85e724993e5fe3c2d30663bd2..435e3200d0aa9d03dfaa21cacd4ff3e65aad712c 100644
--- a/src/space.h
+++ b/src/space.h
@@ -72,9 +72,6 @@ struct space {
/*! Are we doing gravity? */
int gravity;
- /*! Total mass in the system */
- double total_mass;
-
/*! Width of the top-level cells. */
double width[3];