diff --git a/src/engine.c b/src/engine.c index 72136782949c51cea111faf9752b08d539cdcc48..ae8a74d5d040579ce08ca5f6a16f587471a06eb7 100644 --- a/src/engine.c +++ b/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."); diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h index 3f97070fbb80c8920a929e6df8d93def6ac6041a..048bc3ee832d4cf9f939c41a0056ca490729f9fa 100644 --- a/src/gravity/Default/gravity.h +++ b/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 } /** diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h index f484b13663059fa5f4f822aa78748fe4ef9d5926..613a958ad2247cdfb12dfda7668db1868e2de8c6 100644 --- a/src/gravity/Default/gravity_part.h +++ b/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; diff --git a/src/runner.c b/src/runner.c index 19c01221d1e261a4a332061350c006f49f3fed79..491ad3269140549f2eb35356e2f3d3baaaab9436 100644 --- a/src/runner.c +++ b/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); } } diff --git a/src/space.h b/src/space.h index 752974c69055ff7c077fa4e38c54cfa48744d249..1f76e7439e65bf8daf1c77bd734c6392b08e13fe 100644 --- a/src/space.h +++ b/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];