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

Restore the check to verify that the multipoles have been correctly construted and received.

parent 3e365581
......@@ -1419,25 +1419,64 @@ void cell_make_multipoles(struct cell *c, integertime_t ti_current) {
c->ti_old_multipole = ti_current;
}
/**
* @brief Recursively verify that the multipoles are the sum of their progenies.
*
* This function does not check whether the multipoles match the particle
* content as we may not have received the particles.
*
* @param c The #cell to recursively search and verify.
*/
void cell_check_foreign_multipole(const struct cell *c) {
#ifdef SWIFT_DEBUG_CHECKS
if (c->split) {
double M_000 = 0.;
long long num_gpart = 0;
for (int k = 0; k < 8; k++) {
const struct cell *cp = c->progeny[k];
if (cp != NULL) {
/* Check the mass */
M_000 += cp->multipole->m_pole.M_000;
/* Check the number of particles */
num_gpart += cp->multipole->m_pole.num_gpart;
/* Now recurse */
cell_check_foreign_multipole(cp);
}
}
if (num_gpart != c->multipole->m_pole.num_gpart)
error("Sum of particles in progenies does not match");
}
#else
error("Calling debugging code without debugging flag activated.");
#endif
}
/**
* @brief Computes the multi-pole brutally and compare to the
* recursively computed one.
*
* @param c Cell to act upon
* @param data Unused parameter
*/
void cell_check_multipole(struct cell *c, void *data) {
void cell_check_multipole(struct cell *c) {
#ifdef SWIFT_DEBUG_CHECKS
struct gravity_tensors ma;
const double tolerance = 1e-3; /* Relative */
return;
/* First recurse */
if (c->split)
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) cell_check_multipole(c->progeny[k], NULL);
if (c->progeny[k] != NULL) cell_check_multipole(c->progeny[k]);
if (c->gcount > 0) {
......@@ -1456,7 +1495,7 @@ void cell_check_multipole(struct cell *c, void *data) {
}
/* Check that the upper limit of r_max is good enough */
if (!(c->multipole->r_max >= ma.r_max)) {
if (!(1.1 * c->multipole->r_max >= ma.r_max)) {
error("Upper-limit r_max=%e too small. Should be >=%e.",
c->multipole->r_max, ma.r_max);
} else if (c->multipole->r_max * c->multipole->r_max >
......
......@@ -515,7 +515,8 @@ int cell_link_gparts(struct cell *c, struct gpart *gparts);
int cell_link_sparts(struct cell *c, struct spart *sparts);
void cell_clean_links(struct cell *c, void *data);
void cell_make_multipoles(struct cell *c, integertime_t ti_current);
void cell_check_multipole(struct cell *c, void *data);
void cell_check_multipole(struct cell *c);
void cell_check_foreign_multipole(const struct cell *c);
void cell_clean(struct cell *c);
void cell_check_part_drift_point(struct cell *c, void *data);
void cell_check_gpart_drift_point(struct cell *c, void *data);
......
......@@ -3969,6 +3969,18 @@ void engine_rebuild(struct engine *e, int clean_smoothing_length_values) {
engine_exchange_cells(e);
#endif
#ifdef SWIFT_DEBUG_CHECKS
long long counter = 0;
/* Let's check that what we received makes sense */
for (int i = 0; i < e->s->nr_cells; ++i) {
const struct gravity_tensors *m = &e->s->multipoles_top[i];
counter += m->m_pole.num_gpart;
}
if (counter != e->total_nr_gparts)
error("Total particles in multipoles inconsistent with engine");
#endif
/* Re-build the tasks. */
engine_maketasks(e);
......@@ -3981,6 +3993,9 @@ void engine_rebuild(struct engine *e, int clean_smoothing_length_values) {
* previously been active on this rank. */
space_check_drift_point(e->s, e->ti_current,
e->policy & engine_policy_self_gravity);
for (int k = 0; k < e->s->nr_local_cells; k++)
cell_check_foreign_multipole(&e->s->cells_top[e->s->local_cells_top[k]]);
#endif
/* Run through the tasks and mark as skip or not. */
......
......@@ -960,7 +960,7 @@ void space_rebuild(struct space *s, int verbose) {
/* Check that the multipole construction went OK */
if (s->gravity)
for (int k = 0; k < s->nr_cells; k++)
cell_check_multipole(&s->cells_top[k], NULL);
cell_check_multipole(&s->cells_top[k]);
#endif
/* Clean up any stray sort indices in the cell buffer. */
......@@ -2311,6 +2311,7 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
if (cells[j]->sort[k] != NULL) free(cells[j]->sort[k]);
struct gravity_tensors *temp = cells[j]->multipole;
bzero(cells[j], sizeof(struct cell));
bzero(temp, sizeof(struct gravity_tensors));
cells[j]->multipole = temp;
cells[j]->nodeID = -1;
if (lock_init(&cells[j]->lock) != 0 || lock_init(&cells[j]->glock) != 0 ||
......
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