Commit e8bec0ca authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Check that the top-level multipole construction went well

parent 14304352
......@@ -2993,6 +2993,18 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
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 */
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.mass;
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);
}
#endif
/* Now time to get ready for the first time-step */
......@@ -3080,6 +3092,20 @@ void engine_step(struct engine *e) {
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that we have the correct total mass in the top-level multipoles */
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.mass;
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);
}
#endif
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
/* Run the brute-force gravity calculation for some gparts */
gravity_exact_force_compute(e->s, e);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment