Compute the total mass in the system and add a variable to the gpart to record...

Compute the total mass in the system and add a variable to the gpart to record how much mass they interacted with via self-gravity.

src/engine.c

@@ -2950,6 +2950,19 @@ 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
+  message("Total mass in the system: %e", e->s->total_mass);
+#endif
+
   /* Now time to get ready for the first time-step */
   if (e->nodeID == 0) message("Running initial fake time-step.");
 

src/gravity/Default/gravity.h

@@ -57,6 +57,10 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
   gp->a_grav[0] = 0.f;
   gp->a_grav[1] = 0.f;
   gp->a_grav[2] = 0.f;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  gp->mass_interacted = 0.;
+#endif
 }
 
 /**

src/gravity/Default/gravity_part.h

@@ -52,6 +52,9 @@ struct gpart {
 
 #ifdef SWIFT_DEBUG_CHECKS
 
+  /* Total mass this gpart interacted with */
+  double mass_interacted;
+
   /* Time of the last drift */
   integertime_t ti_drift;
 

src/runner.c

@@ -1332,13 +1332,7 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
 
       /* Get a handle on the part. */
       struct part *restrict p = &parts[k];
-
-      if (part_is_active(p, e)) {
-
-        /* First, finish the force loop */
-        hydro_end_force(p);
-        if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
-      }
+      if (part_is_active(p, e)) hydro_end_force(p);
     }
 
     /* Loop over the g-particles in this cell. */
@@ -1346,24 +1340,26 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
 
       /* Get a handle on the gpart. */
       struct gpart *restrict gp = &gparts[k];
+      if (gpart_is_active(gp, e)) gravity_end_force(gp, const_G);
 
-      if (gp->type == swift_type_dark_matter) {
-        if (gpart_is_active(gp, e)) gravity_end_force(gp, const_G);
+#ifdef SWIFT_DEBUG_CHECKS
+      if (e->policy & engine_policy_self_gravity) {
+        gp->mass_interacted += gp->mass;
+        if (gp->mass_interacted != e->s->total_mass)
+          error(
+              "g-particle did not interact gravitationally with all other "
+              "particles gp->mass_interacted=%e, total_mass=%e",
+              gp->mass_interacted, e->s->total_mass);
       }
+#endif
     }
 
     /* Loop over the star particles in this cell. */
     for (int k = 0; k < scount; k++) {
 
-      /* Get a handle on the part. */
+      /* Get a handle on the spart. */
       struct spart *restrict sp = &sparts[k];
-
-      if (spart_is_active(sp, e)) {
-
-        /* First, finish the force loop */
-        star_end_force(sp);
-        gravity_end_force(sp->gpart, const_G);
-      }
+      if (spart_is_active(sp, e)) star_end_force(sp);
     }
   }
 

src/space.h

@@ -72,6 +72,9 @@ 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];