Commit 56d1a511 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'ic-checks' into 'master'

Add checks for particles that are at same position and updates to the h_max of top-level cells.

Investigating #324 has turned up a couple of issues that we probably should
handle better. The first is particles with the same position. These are not 
allowed. The second is a mismatch between the cells h_max values and those
of the particles when running the ghost tasks for the first time and when
the initial estimates of h are small, so we get too many top-level cells 
initially and a regrid is needed (if this isn't done then the particles
test to be moved outside their cell and we need to rebuild, even though
a rebuild has just been done).

This fixes both by reporting an error when particles are found in the
same location, and updates the top-level cells h_max to those of the
particles (note this must be followed by a regrid as any sub-cells
will have the wrong h_max).

Also now reports if gparts are at the same location. That is inefficient
not a show stopper.

See merge request !370
parents 7bf76a88 eef6e5de
......@@ -3360,6 +3360,67 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
/* Recover the (integer) end of the next time-step */
engine_collect_timestep_and_rebuild(e, 1);
/* Check if any particles have the same position. This is not
* allowed (/0) so we abort.*/
if (s->nr_parts > 0) {
/* Sorting should put the same positions next to each other... */
int failed = 0;
double *prev_x = s->parts[0].x;
for (size_t k = 1; k < s->nr_parts; k++) {
if (prev_x[0] == s->parts[k].x[0] &&
prev_x[1] == s->parts[k].x[1] &&
prev_x[2] == s->parts[k].x[2]) {
if (e->verbose)
message("Two particles occupy location: %f %f %f",
prev_x[0], prev_x[1], prev_x[2]);
prev_x = s->parts[k].x;
if (failed > 0)
error("Have %d particle pairs with the same locations.\n"
"Cannot continue", failed);
/* Also check any gparts. This is not supposed to be fatal so only warn. */
if (s->nr_gparts > 0) {
int failed = 0;
double *prev_x = s->gparts[0].x;
for (size_t k = 1; k < s->nr_gparts; k++) {
if (prev_x[0] == s->gparts[k].x[0] &&
prev_x[1] == s->gparts[k].x[1] &&
prev_x[2] == s->gparts[k].x[2]) {
if (e->verbose)
message("Two gparts occupy location: %f %f %f / %f %f %f",
prev_x[0], prev_x[1], prev_x[2],
s->gparts[k].x[0], s->gparts[k].x[1], s->gparts[k].x[2] );
prev_x = s->gparts[k].x;
if (failed > 0)
message("WARNING: found %d gpart pairs at the same location. "
"That is not optimal", failed);
/* Check the top-level cell h_max matches the particles as these can be
* updated in the the ghost tasks (only a problem if the ICs estimates for h
* are too small). Note this must be followed by a rebuild as sub-cells will
* not be updated until that is done. */
if (s->cells_top != NULL && s->nr_parts > 0) {
for (int i = 0; i < s->nr_cells; i++) {
struct cell *c = &s->cells_top[i];
if (c->nodeID == engine_rank) {
float part_h_max = c->parts[0].h;
for (int k = 1; k < c->count; k++) {
if (c->parts[k].h > part_h_max) part_h_max = c->parts[k].h;
c->h_max = max(part_h_max, c->h_max);
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