Skip to content

SWIFTsim
SWIFTsim
Commit 86db3dab authored by Peter W. Draper's avatar Peter W. Draper
Browse files
Options

Merge branch 'faster_rebuilds' into 'master' 

Updated cell splitting strategy

See merge request !493
parents 2643736b 7cb536a4
Hide whitespace changes
Inline Side-by-side
Showing with 337 additions and 623 deletions
src/cell.c
View file @ 86db3dab
......@@ -835,11 +835,15 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
memswap(&parts[j], &part, sizeof(struct part));
memswap(&xparts[j], &xpart, sizeof(struct xpart));
memswap(&buff[j], &temp_buff, sizeof(struct cell_buff));
if (parts[j].gpart)
parts[j].gpart->id_or_neg_offset = -(j + parts_offset);
bid = temp_buff.ind;
}
parts[k] = part;
xparts[k] = xpart;
buff[k] = temp_buff;
if (parts[k].gpart)
parts[k].gpart->id_or_neg_offset = -(k + parts_offset);
}
bucket_count[bid]++;
}
......@@ -852,10 +856,6 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
c->progeny[k]->xparts = &c->xparts[bucket_offset[k]];
}
/* Re-link the gparts. */
if (count > 0 && gcount > 0)
part_relink_gparts_to_parts(parts, count, parts_offset);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that the buffs are OK. */
for (int k = 1; k < count; k++) {
......@@ -952,10 +952,14 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
}
memswap(&sparts[j], &spart, sizeof(struct spart));
memswap(&sbuff[j], &temp_buff, sizeof(struct cell_buff));
if (sparts[j].gpart)
sparts[j].gpart->id_or_neg_offset = -(j + sparts_offset);
bid = temp_buff.ind;
}
sparts[k] = spart;
sbuff[k] = temp_buff;
if (sparts[k].gpart)
sparts[k].gpart->id_or_neg_offset = -(k + sparts_offset);
}
bucket_count[bid]++;
}
......@@ -967,10 +971,6 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
c->progeny[k]->sparts = &c->sparts[bucket_offset[k]];
}
/* Re-link the gparts. */
if (scount > 0 && gcount > 0)
part_relink_gparts_to_sparts(sparts, scount, sparts_offset);
/* Finally, do the same song and dance for the gparts. */
for (int k = 0; k < 8; k++) bucket_count[k] = 0;
......@@ -1005,10 +1005,23 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
}
memswap(&gparts[j], &gpart, sizeof(struct gpart));
memswap(&gbuff[j], &temp_buff, sizeof(struct cell_buff));
if (gparts[j].type == swift_type_gas) {
parts[-gparts[j].id_or_neg_offset - parts_offset].gpart =
&gparts[j];
} else if (gparts[j].type == swift_type_star) {
sparts[-gparts[j].id_or_neg_offset - sparts_offset].gpart =
&gparts[j];
}
bid = temp_buff.ind;
}
gparts[k] = gpart;
gbuff[k] = temp_buff;
if (gparts[k].type == swift_type_gas) {
parts[-gparts[k].id_or_neg_offset - parts_offset].gpart = &gparts[k];
} else if (gparts[k].type == swift_type_star) {
sparts[-gparts[k].id_or_neg_offset - sparts_offset].gpart =
&gparts[k];
}
}
bucket_count[bid]++;
}
......@@ -1019,14 +1032,6 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
c->progeny[k]->gcount = bucket_count[k];
c->progeny[k]->gparts = &c->gparts[bucket_offset[k]];
}
/* Re-link the parts. */
if (count > 0 && gcount > 0)
part_relink_parts_to_gparts(gparts, gcount, parts - parts_offset);
/* Re-link the sparts. */
if (scount > 0 && gcount > 0)
part_relink_sparts_to_gparts(gparts, gcount, sparts - sparts_offset);
}
/**
......
src/engine.c
View file @ 86db3dab
......@@ -565,7 +565,8 @@ void engine_redistribute(struct engine *e) {
/* Sort the particles according to their cell index. */
if (s->nr_parts > 0)
space_parts_sort(s, dest, s->nr_parts, 0, nr_nodes - 1, e->verbose);
space_parts_sort(s->parts, s->xparts, dest, &counts[nodeID * nr_nodes],
nr_nodes, 0);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the part have been sorted correctly. */
......@@ -656,7 +657,8 @@ void engine_redistribute(struct engine *e) {
/* Sort the particles according to their cell index. */
if (s->nr_sparts > 0)
space_sparts_sort(s, s_dest, s->nr_sparts, 0, nr_nodes - 1, e->verbose);
space_sparts_sort(s->sparts, s_dest, &s_counts[nodeID * nr_nodes], nr_nodes,
0);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the spart have been sorted correctly. */
......@@ -748,7 +750,8 @@ void engine_redistribute(struct engine *e) {
/* Sort the gparticles according to their cell index. */
if (s->nr_gparts > 0)
space_gparts_sort(s, g_dest, s->nr_gparts, 0, nr_nodes - 1, e->verbose);
space_gparts_sort(s->gparts, s->parts, s->sparts, g_dest,
&g_counts[nodeID * nr_nodes], nr_nodes);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the gpart have been sorted correctly. */
......@@ -3678,6 +3681,11 @@ void engine_rebuild(struct engine *e, int clean_smoothing_length_values) {
/* Re-build the space. */
space_rebuild(e->s, e->verbose);
#ifdef SWIFT_DEBUG_CHECKS
part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->nr_parts,
e->s->nr_gparts, e->s->nr_sparts, e->verbose);
#endif
/* Initial cleaning up session ? */
if (clean_smoothing_length_values) space_sanitize(e->s);
......@@ -4151,7 +4159,7 @@ void engine_first_init_particles(struct engine *e) {
const ticks tic = getticks();
/* Set the particles in a state where they are ready for a run */
/* Set the particles in a state where they are ready for a run. */
space_first_init_parts(e->s, e->chemistry, e->cooling_func);
space_first_init_gparts(e->s, e->gravity_properties);
space_first_init_sparts(e->s);
......
src/part.c
View file @ 86db3dab
......@@ -94,6 +94,26 @@ void part_relink_sparts_to_gparts(struct gpart *gparts, size_t N,
}
}
/**
* @brief Re-link both the #part%s and #spart%s associated with the list of
* #gpart%s.
*
* @param gparts The list of #gpart.
* @param N The number of particles to re-link;
* @param parts The global #part array in which to find the #gpart offsets.
* @param sparts The global #spart array in which to find the #gpart offsets.
*/
void part_relink_all_parts_to_gparts(struct gpart *gparts, size_t N,
struct part *parts, struct spart *sparts) {
for (size_t k = 0; k < N; k++) {
if (gparts[k].type == swift_type_gas) {
parts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
} else if (gparts[k].type == swift_type_star) {
sparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
}
}
}
/**
* @brief Verifies that the #gpart, #part and #spart are correctly linked
* together
......@@ -128,19 +148,19 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
/* Check that it is linked */
if (gparts[k].id_or_neg_offset > 0)
error("Gas gpart not linked to anything !");
error("Gas gpart not linked to anything!");
/* Find its link */
const struct part *part = &parts[-gparts[k].id_or_neg_offset];
/* Check the reverse link */
if (part->gpart != &gparts[k]) error("Linking problem !");
if (part->gpart != &gparts[k]) error("Linking problem!");
/* Check that the particles are at the same place */
if (gparts[k].x[0] != part->x[0] || gparts[k].x[1] != part->x[1] ||
gparts[k].x[2] != part->x[2])
error(
"Linked particles are not at the same position !\n"
"Linked particles are not at the same position!\n"
"gp->x=[%e %e %e] p->x=[%e %e %e] diff=[%e %e %e]",
gparts[k].x[0], gparts[k].x[1], gparts[k].x[2], part->x[0],
part->x[1], part->x[2], gparts[k].x[0] - part->x[0],
......
src/part.h
View file @ 86db3dab
......@@ -80,6 +80,8 @@ void part_relink_parts_to_gparts(struct gpart *gparts, size_t N,
struct part *parts);
void part_relink_sparts_to_gparts(struct gpart *gparts, size_t N,
struct spart *sparts);
void part_relink_all_parts_to_gparts(struct gpart *gparts, size_t N,
struct part *parts, struct spart *sparts);
void part_verify_links(struct part *parts, struct gpart *gparts,
struct spart *sparts, size_t nr_parts, size_t nr_gparts,
size_t nr_sparts, int verbose);
......
src/space.c
View file @ 86db3dab
......@@ -101,6 +101,7 @@ struct index_data {
struct space *s;
struct cell *cells;
int *ind;
int *cell_counts;
};
/**
......@@ -580,26 +581,33 @@ void space_rebuild(struct space *s, int verbose) {
an index that is larger than the number of particles to avoid
re-allocating after shuffling. */
const size_t ind_size = s->size_parts + 100;
int *ind;
if ((ind = (int *)malloc(sizeof(int) * ind_size)) == NULL)
error("Failed to allocate temporary particle indices.");
if (s->size_parts > 0) space_parts_get_cell_index(s, ind, cells_top, verbose);
int *ind = (int *)malloc(sizeof(int) * ind_size);
if (ind == NULL) error("Failed to allocate temporary particle indices.");
int *cell_part_counts = (int *)calloc(sizeof(int), s->nr_cells);
if (cell_part_counts == NULL)
error("Failed to allocate cell part count buffer.");
if (s->size_parts > 0)
space_parts_get_cell_index(s, ind, cell_part_counts, cells_top, verbose);
/* Run through the gravity particles and get their cell index. */
const size_t gind_size = s->size_gparts + 100;
int *gind;
if ((gind = (int *)malloc(sizeof(int) * gind_size)) == NULL)
error("Failed to allocate temporary g-particle indices.");
int *gind = (int *)malloc(sizeof(int) * gind_size);
if (gind == NULL) error("Failed to allocate temporary g-particle indices.");
int *cell_gpart_counts = (int *)calloc(sizeof(int), s->nr_cells);
if (cell_gpart_counts == NULL)
error("Failed to allocate cell gpart count buffer.");
if (s->size_gparts > 0)
space_gparts_get_cell_index(s, gind, cells_top, verbose);
space_gparts_get_cell_index(s, gind, cell_gpart_counts, cells_top, verbose);
/* Run through the star particles and get their cell index. */
const size_t sind_size = s->size_sparts + 100;
int *sind;
if ((sind = (int *)malloc(sizeof(int) * sind_size)) == NULL)
error("Failed to allocate temporary s-particle indices.");
int *sind = (int *)malloc(sizeof(int) * sind_size);
if (sind == NULL) error("Failed to allocate temporary s-particle indices.");
int *cell_spart_counts = (int *)calloc(sizeof(int), s->nr_cells);
if (cell_spart_counts == NULL)
error("Failed to allocate cell gpart count buffer.");
if (s->size_sparts > 0)
space_sparts_get_cell_index(s, sind, cells_top, verbose);
space_sparts_get_cell_index(s, sind, cell_spart_counts, cells_top, verbose);
#ifdef WITH_MPI
const int local_nodeID = s->e->nodeID;
......@@ -609,9 +617,7 @@ void space_rebuild(struct space *s, int verbose) {
if (cells_top[ind[k]].nodeID != local_nodeID) {
nr_parts -= 1;
/* Swap the particle */
const struct part tp = s->parts[k];
s->parts[k] = s->parts[nr_parts];
s->parts[nr_parts] = tp;
memswap(&s->parts[k], &s->parts[nr_parts], sizeof(struct part));
/* Swap the link with the gpart */
if (s->parts[k].gpart != NULL) {
s->parts[k].gpart->id_or_neg_offset = -k;
......@@ -620,13 +626,9 @@ void space_rebuild(struct space *s, int verbose) {
s->parts[nr_parts].gpart->id_or_neg_offset = -nr_parts;
}
/* Swap the xpart */
const struct xpart txp = s->xparts[k];
s->xparts[k] = s->xparts[nr_parts];
s->xparts[nr_parts] = txp;
memswap(&s->xparts[k], &s->xparts[nr_parts], sizeof(struct xpart));
/* Swap the index */
const int t = ind[k];
ind[k] = ind[nr_parts];
ind[nr_parts] = t;
memswap(&ind[k], &ind[nr_parts], sizeof(int));
} else {
/* Increment when not exchanging otherwise we need to retest "k".*/
k++;
......@@ -652,9 +654,7 @@ void space_rebuild(struct space *s, int verbose) {
if (cells_top[sind[k]].nodeID != local_nodeID) {
nr_sparts -= 1;
/* Swap the particle */
const struct spart tp = s->sparts[k];
s->sparts[k] = s->sparts[nr_sparts];
s->sparts[nr_sparts] = tp;
memswap(&s->sparts[k], &s->sparts[nr_sparts], sizeof(struct spart));
/* Swap the link with the gpart */
if (s->sparts[k].gpart != NULL) {
s->sparts[k].gpart->id_or_neg_offset = -k;
......@@ -663,9 +663,7 @@ void space_rebuild(struct space *s, int verbose) {
s->sparts[nr_sparts].gpart->id_or_neg_offset = -nr_sparts;
}
/* Swap the index */
const int t = sind[k];
sind[k] = sind[nr_sparts];
sind[nr_sparts] = t;
memswap(&sind[k], &sind[nr_sparts], sizeof(int));
} else {
/* Increment when not exchanging otherwise we need to retest "k".*/
k++;
......@@ -691,9 +689,7 @@ void space_rebuild(struct space *s, int verbose) {
if (cells_top[gind[k]].nodeID != local_nodeID) {
nr_gparts -= 1;
/* Swap the particle */
const struct gpart tp = s->gparts[k];
s->gparts[k] = s->gparts[nr_gparts];
s->gparts[nr_gparts] = tp;
memswap(&s->gparts[k], &s->gparts[nr_gparts], sizeof(struct gpart));
/* Swap the link with part/spart */
if (s->gparts[k].type == swift_type_gas) {
s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
......@@ -708,9 +704,7 @@ void space_rebuild(struct space *s, int verbose) {
&s->gparts[nr_gparts];
}
/* Swap the index */
const int t = gind[k];
gind[k] = gind[nr_gparts];
gind[nr_gparts] = t;
memswap(&gind[k], &gind[nr_gparts], sizeof(int));
} else {
/* Increment when not exchanging otherwise we need to retest "k".*/
k++;
......@@ -745,6 +739,15 @@ void space_rebuild(struct space *s, int verbose) {
s->nr_gparts = nr_gparts + nr_gparts_exchanged;
s->nr_sparts = nr_sparts + nr_sparts_exchanged;
/* Clear non-local cell counts. */
for (int k = 0; k < s->nr_cells; k++) {
if (s->cells_top[k].nodeID != local_nodeID) {
cell_part_counts[k] = 0;
cell_spart_counts[k] = 0;
cell_gpart_counts[k] = 0;
}
}
/* Re-allocate the index array for the parts if needed.. */
if (s->nr_parts + 1 > ind_size) {
int *ind_new;
......@@ -773,6 +776,7 @@ void space_rebuild(struct space *s, int verbose) {
const struct part *const p = &s->parts[k];
ind[k] =
cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
cell_part_counts[ind[k]]++;
#ifdef SWIFT_DEBUG_CHECKS
if (cells_top[ind[k]].nodeID != local_nodeID)
error("Received part that does not belong to me (nodeID=%i).",
......@@ -786,6 +790,7 @@ void space_rebuild(struct space *s, int verbose) {
const struct spart *const sp = &s->sparts[k];
sind[k] =
cell_getid(cdim, sp->x[0] * ih[0], sp->x[1] * ih[1], sp->x[2] * ih[2]);
cell_spart_counts[sind[k]]++;
#ifdef SWIFT_DEBUG_CHECKS
if (cells_top[sind[k]].nodeID != local_nodeID)
error("Received s-part that does not belong to me (nodeID=%i).",
......@@ -794,11 +799,12 @@ void space_rebuild(struct space *s, int verbose) {
}
nr_sparts = s->nr_sparts;
#endif /* WITH_MPI */
#endif // WITH_MPI
/* Sort the parts according to their cells. */
if (nr_parts > 0)
space_parts_sort(s, ind, nr_parts, 0, s->nr_cells - 1, verbose);
space_parts_sort(s->parts, s->xparts, ind, cell_part_counts, s->nr_cells,
0);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the part have been sorted correctly. */
......@@ -825,7 +831,7 @@ void space_rebuild(struct space *s, int verbose) {
/* Sort the sparts according to their cells. */
if (nr_sparts > 0)
space_sparts_sort(s, sind, nr_sparts, 0, s->nr_cells - 1, verbose);
space_sparts_sort(s->sparts, sind, cell_spart_counts, s->nr_cells, 0);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the spart have been sorted correctly. */
......@@ -850,12 +856,6 @@ void space_rebuild(struct space *s, int verbose) {
}
#endif
/* Re-link the gparts to their (s-)particles. */
if (nr_parts > 0 && nr_gparts > 0)
part_relink_gparts_to_parts(s->parts, nr_parts, 0);
if (nr_sparts > 0 && nr_gparts > 0)
part_relink_gparts_to_sparts(s->sparts, nr_sparts, 0);
/* Extract the cell counts from the sorted indices. */
size_t last_index = 0;
ind[nr_parts] = s->nr_cells; // sentinel.
......@@ -878,7 +878,9 @@ void space_rebuild(struct space *s, int verbose) {
/* We no longer need the indices as of here. */
free(ind);
free(cell_part_counts);
free(sind);
free(cell_spart_counts);
#ifdef WITH_MPI
......@@ -897,7 +899,7 @@ void space_rebuild(struct space *s, int verbose) {
const struct gpart *const p = &s->gparts[k];
gind[k] =
cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
cell_gpart_counts[gind[k]]++;
#ifdef SWIFT_DEBUG_CHECKS
if (cells_top[gind[k]].nodeID != s->e->nodeID)
error("Received g-part that does not belong to me (nodeID=%i).",
......@@ -906,11 +908,12 @@ void space_rebuild(struct space *s, int verbose) {
}
nr_gparts = s->nr_gparts;
#endif /* WITH_MPI */
#endif // WITH_MPI
/* Sort the gparts according to their cells. */
if (nr_gparts > 0)
space_gparts_sort(s, gind, nr_gparts, 0, s->nr_cells - 1, verbose);
space_gparts_sort(s->gparts, s->parts, s->sparts, gind, cell_gpart_counts,
s->nr_cells);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the gpart have been sorted correctly. */
......@@ -935,14 +938,6 @@ void space_rebuild(struct space *s, int verbose) {
}
#endif
/* Re-link the parts. */
if (nr_parts > 0 && nr_gparts > 0)
part_relink_parts_to_gparts(s->gparts, nr_gparts, s->parts);
/* Re-link the sparts. */
if (nr_sparts > 0 && nr_gparts > 0)
part_relink_sparts_to_gparts(s->gparts, nr_gparts, s->sparts);
/* Extract the cell counts from the sorted indices. */
size_t last_gindex = 0;
gind[nr_gparts] = s->nr_cells;
......@@ -955,6 +950,7 @@ void space_rebuild(struct space *s, int verbose) {
/* We no longer need the indices as of here. */
free(gind);
free(cell_gpart_counts);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the links are correct */
......@@ -999,6 +995,9 @@ void space_rebuild(struct space *s, int verbose) {
cell_check_multipole(&s->cells_top[k], NULL);
#endif
/* Clean up any stray sort indices in the cell buffer. */
space_free_buff_sort_indices(s);
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
......@@ -1082,6 +1081,12 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts,
const double ih_y = s->iwidth[1];
const double ih_z = s->iwidth[2];
/* Init the local count buffer. */
int *cell_counts = (int *)calloc(sizeof(int), s->nr_cells);
if (cell_counts == NULL)
error("Failed to allocate temporary cell count buffer.");
/* Loop over the parts. */
for (int k = 0; k < nr_parts; k++) {
/* Get the particle */
......@@ -1100,6 +1105,7 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts,
const int index =
cell_getid(cdim, pos_x * ih_x, pos_y * ih_y, pos_z * ih_z);
ind[k] = index;
cell_counts[index]++;
#ifdef SWIFT_DEBUG_CHECKS
if (index < 0 || index >= cdim[0] * cdim[1] * cdim[2])
......@@ -1117,6 +1123,11 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts,
p->x[1] = pos_y;
p->x[2] = pos_z;
}
/* Write the counts back to the global array. */
for (int k = 0; k < s->nr_cells; k++)
if (cell_counts[k]) atomic_add(&data->cell_counts[k], cell_counts[k]);
free(cell_counts);
}
/**
......@@ -1144,6 +1155,11 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts,
const double ih_y = s->iwidth[1];
const double ih_z = s->iwidth[2];
/* Init the local count buffer. */
int *cell_counts = (int *)calloc(sizeof(int), s->nr_cells);
if (cell_counts == NULL)
error("Failed to allocate temporary cell count buffer.");
for (int k = 0; k < nr_gparts; k++) {
/* Get the particle */
......@@ -1162,6 +1178,7 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts,
const int index =
cell_getid(cdim, pos_x * ih_x, pos_y * ih_y, pos_z * ih_z);
ind[k] = index;
cell_counts[index]++;
#ifdef SWIFT_DEBUG_CHECKS
if (index < 0 || index >= cdim[0] * cdim[1] * cdim[2])
......@@ -1179,6 +1196,11 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts,
gp->x[1] = pos_y;
gp->x[2] = pos_z;
}
/* Write the counts back to the global array. */
for (int k = 0; k < s->nr_cells; k++)
if (cell_counts[k]) atomic_add(&data->cell_counts[k], cell_counts[k]);
free(cell_counts);
}
/**
......@@ -1206,6 +1228,11 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
const double ih_y = s->iwidth[1];
const double ih_z = s->iwidth[2];
/* Init the local count buffer. */
int *cell_counts = (int *)calloc(sizeof(int), s->nr_cells);
if (cell_counts == NULL)
error("Failed to allocate temporary cell count buffer.");
for (int k = 0; k < nr_sparts; k++) {
/* Get the particle */
......@@ -1224,6 +1251,7 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
const int index =
cell_getid(cdim, pos_x * ih_x, pos_y * ih_y, pos_z * ih_z);
ind[k] = index;
cell_counts[index]++;
#ifdef SWIFT_DEBUG_CHECKS
if (index < 0 || index >= cdim[0] * cdim[1] * cdim[2])
......@@ -1241,6 +1269,11 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
sp->x[1] = pos_y;
sp->x[2] = pos_z;
}
/* Write the counts back to the global array. */
for (int k = 0; k < s->nr_cells; k++)
if (cell_counts[k]) atomic_add(&data->cell_counts[k], cell_counts[k]);
free(cell_counts);
}
/**
......@@ -1248,11 +1281,12 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
*
* @param s The #space.
* @param ind The array of indices to fill.
* @param cell_counts The cell counters to update.
* @param cells The array of #cell to update.
* @param verbose Are we talkative ?
*/
void space_parts_get_cell_index(struct space *s, int *ind, struct cell *cells,
int verbose) {
void space_parts_get_cell_index(struct space *s, int *ind, int *cell_counts,
struct cell *cells, int verbose) {
const ticks tic = getticks();
......@@ -1261,6 +1295,7 @@ void space_parts_get_cell_index(struct space *s, int *ind, struct cell *cells,
data.s = s;
data.cells = cells;
data.ind = ind;
data.cell_counts = cell_counts;
threadpool_map(&s->e->threadpool, space_parts_get_cell_index_mapper, s->parts,
s->nr_parts, sizeof(struct part), 0, &data);
......@@ -1275,11 +1310,12 @@ void space_parts_get_cell_index(struct space *s, int *ind, struct cell *cells,
*
* @param s The #space.
* @param gind The array of indices to fill.
* @param cell_counts The cell counters to update.
* @param cells The array of #cell to update.
* @param verbose Are we talkative ?
*/
void space_gparts_get_cell_index(struct space *s, int *gind, struct cell *cells,
int verbose) {
void space_gparts_get_cell_index(struct space *s, int *gind, int *cell_counts,
struct cell *cells, int verbose) {
const ticks tic = getticks();
......@@ -1288,6 +1324,7 @@ void space_gparts_get_cell_index(struct space *s, int *gind, struct cell *cells,
data.s = s;
data.cells = cells;
data.ind = gind;
data.cell_counts = cell_counts;
threadpool_map(&s->e->threadpool, space_gparts_get_cell_index_mapper,
s->gparts, s->nr_gparts, sizeof(struct gpart), 0, &data);
......@@ -1302,11 +1339,12 @@ void space_gparts_get_cell_index(struct space *s, int *gind, struct cell *cells,
*
* @param s The #space.
* @param sind The array of indices to fill.
* @param cell_counts The cell counters to update.
* @param cells The array of #cell to update.
* @param verbose Are we talkative ?
*/
void space_sparts_get_cell_index(struct space *s, int *sind, struct cell *cells,
int verbose) {
void space_sparts_get_cell_index(struct space *s, int *sind, int *cell_counts,
struct cell *cells, int verbose) {
const ticks tic = getticks();
......@@ -1315,6 +1353,7 @@ void space_sparts_get_cell_index(struct space *s, int *sind, struct cell *cells,
data.s = s;
data.cells = cells;
data.ind = sind;
data.cell_counts = cell_counts;
threadpool_map(&s->e->threadpool, space_sparts_get_cell_index_mapper,
s->sparts, s->nr_sparts, sizeof(struct spart), 0, &data);
......@@ -1328,551 +1367,176 @@ void space_sparts_get_cell_index(struct space *s, int *sind, struct cell *cells,
* @brief Sort the particles and condensed particles according to the given
* indices.
*
* @param s The #space.
* @param parts The array of #part to sort.
* @param xparts The corresponding #xpart array to sort as well.
* @param ind The indices with respect to which the parts are sorted.
* @param N The number of parts
* @param min Lowest index.
* @param max highest index.
* @param verbose Are we talkative ?
* @param counts Number of particles per index.
* @param num_bins Total number of bins (length of count).
* @param parts_offset Offset of the #part array from the global #part array.
*/
void space_parts_sort(struct space *s, int *ind, size_t N, int min, int max,
int verbose) {
const ticks tic = getticks();
/* Populate a parallel_sort structure with the input data */
struct parallel_sort sort_struct;
sort_struct.parts = s->parts;
sort_struct.xparts = s->xparts;
sort_struct.ind = ind;
sort_struct.stack_size = 2 * (max - min + 1) + 10 + s->e->nr_threads;
if ((sort_struct.stack = (struct qstack *)malloc(
sizeof(struct qstack) * sort_struct.stack_size)) == NULL)
error("Failed to allocate sorting stack.");
for (unsigned int i = 0; i < sort_struct.stack_size; i++)
sort_struct.stack[i].ready = 0;
/* Add the first interval. */
sort_struct.stack[0].i = 0;
sort_struct.stack[0].j = N - 1;
sort_struct.stack[0].min = min;
sort_struct.stack[0].max = max;
sort_struct.stack[0].ready = 1;
sort_struct.first = 0;
sort_struct.last = 1;
sort_struct.waiting = 1;
/* Launch the sorting tasks with a stride of zero such that the same
map data is passed to each thread. */
threadpool_map(&s->e->threadpool, space_parts_sort_mapper, &sort_struct,
s->e->threadpool.num_threads, 0, 1, NULL);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify space_sort_struct. */
for (size_t i = 1; i < N; i++)
if (ind[i - 1] > ind[i])
error("Sorting failed (ind[%zu]=%i,ind[%zu]=%i), min=%i, max=%i.", i - 1,
ind[i - 1], i, ind[i], min, max);
if (s->e->nodeID == 0 || verbose) message("Sorting succeeded.");
#endif
/* Clean up. */
free(sort_struct.stack);
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
}
void space_parts_sort_mapper(void *map_data, int num_elements,
void *extra_data) {
/* Unpack the mapping data. */
struct parallel_sort *sort_struct = (struct parallel_sort *)map_data;
/* Pointers to the sorting data. */
int *ind = sort_struct->ind;
struct part *parts = sort_struct->parts;
struct xpart *xparts = sort_struct->xparts;
/* Main loop. */
while (sort_struct->waiting) {
/* Grab an interval off the queue. */
int qid = atomic_inc(&sort_struct->first) % sort_struct->stack_size;
/* Wait for the entry to be ready, or for the sorting do be done. */
while (!sort_struct->stack[qid].ready)
if (!sort_struct->waiting) return;
/* Get the stack entry. */
ptrdiff_t i = sort_struct->stack[qid].i;
ptrdiff_t j = sort_struct->stack[qid].j;
int min = sort_struct->stack[qid].min;
int max = sort_struct->stack[qid].max;
sort_struct->stack[qid].ready = 0;
/* Loop over sub-intervals. */
while (1) {
/* Bring beer. */
const int pivot = (min + max) / 2;
/* message("Working on interval [%i,%i] with min=%i, max=%i, pivot=%i.",
i, j, min, max, pivot); */
/* One pass of QuickSort's partitioning. */
ptrdiff_t ii = i;
ptrdiff_t jj = j;
while (ii < jj) {
while (ii <= j && ind[ii] <= pivot) ii++;
while (jj >= i && ind[jj] > pivot) jj--;
if (ii < jj) {