diff --git a/src/minmax.h b/src/minmax.h index 9d92cd71d849dba615fdb05bc342014e0593d989..a53093663c79cf4280d136747663552e49c7f1b2 100644 --- a/src/minmax.h +++ b/src/minmax.h @@ -43,4 +43,18 @@ _a > _b ? _a : _b; \ }) +/** + * @brief Limits the value of x to be between a and b + * + * Only wraps once. If x > 2b, the returned value will be larger than b. + * Similarly for x < -b. + */ +#define box_wrap(x, a, b) \ + ({ \ + const __typeof__(x) _x = (x); \ + const __typeof__(a) _a = (a); \ + const __typeof__(b) _b = (b); \ + _x < _a ? (_x + _b) : ((_x > _b) ? (_x - _b) : _x); \ + }) + #endif /* SWIFT_MINMAX_H */ diff --git a/src/space.c b/src/space.c index fd9c39dd022b32b1654c0e1023330468363f76dc..b40f5f474b1b39914cd3dda454df52b00c1d5d33 100644 --- a/src/space.c +++ b/src/space.c @@ -738,7 +738,7 @@ void space_sanitize(struct space *s) { * @brief #threadpool mapper function to compute the particle cell indices. * * @param map_data Pointer towards the particles. - * @param num_parts The number of particles to treat. + * @param nr_parts The number of particles to treat. * @param extra_data Pointers to the space and index list */ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts, @@ -748,7 +748,7 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts, struct part *restrict parts = (struct part *)map_data; struct index_data *data = (struct index_data *)extra_data; struct space *s = data->s; - int *ind = data->ind; + int *ind = data->ind + (ptrdiff_t)(parts - s->parts); struct cell *cells = data->cells; /* Get some constants */ @@ -762,18 +762,14 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts, struct part *restrict p = &parts[k]; /* Put it back into the simulation volume */ - for (int j = 0; j < 3; j++) - if (p->x[j] < 0.0) - p->x[j] += dim[j]; - else if (p->x[j] >= dim[j]) - p->x[j] -= dim[j]; + p->x[0] = box_wrap(p->x[0], 0.0, dim[0]); + p->x[1] = box_wrap(p->x[1], 0.0, dim[1]); + p->x[2] = box_wrap(p->x[2], 0.0, dim[2]); /* Get its cell index */ const int index = cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]); - - /* Save the index at the right place */ - *(ind + (ptrdiff_t)(p - s->parts)) = index; + ind[k] = index; /* Tell the cell it has a new member */ atomic_inc(&(cells[index].count)); @@ -784,7 +780,7 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts, * @brief #threadpool mapper function to compute the g-particle cell indices. * * @param map_data Pointer towards the g-particles. - * @param num_parts The number of g-particles to treat. + * @param nr_gparts The number of g-particles to treat. * @param extra_data Pointers to the space and index list */ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts, @@ -794,7 +790,7 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts, struct gpart *restrict gparts = (struct gpart *)map_data; struct index_data *data = (struct index_data *)extra_data; struct space *s = data->s; - int *ind = data->ind; + int *ind = data->ind + (ptrdiff_t)(gparts - s->gparts); struct cell *cells = data->cells; /* Get some constants */ @@ -808,18 +804,14 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts, struct gpart *restrict gp = &gparts[k]; /* Put it back into the simulation volume */ - for (int j = 0; j < 3; j++) - if (gp->x[j] < 0.0) - gp->x[j] += dim[j]; - else if (gp->x[j] >= dim[j]) - gp->x[j] -= dim[j]; + gp->x[0] = box_wrap(gp->x[0], 0.0, dim[0]); + gp->x[1] = box_wrap(gp->x[1], 0.0, dim[1]); + gp->x[2] = box_wrap(gp->x[2], 0.0, dim[2]); /* Get its cell index */ const int index = cell_getid(cdim, gp->x[0] * ih[0], gp->x[1] * ih[1], gp->x[2] * ih[2]); - - /* Save the index at the right place */ - *(ind + (ptrdiff_t)(gp - s->gparts)) = index; + ind[k] = index; /* Tell the cell it has a new member */ atomic_inc(&(cells[index].gcount));