diff --git a/src/Makefile.am b/src/Makefile.am
index 1b9d3b4ae9770b9081ba384851fe6cbf1f2ad688..d273dde02dd186de3697fbc42deb5b49c951941f 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -59,7 +59,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \
units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
- dimension.h equation_of_state.h active.h \
+ dimension.h equation_of_state.h active.h part_type.h \
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
diff --git a/src/cell.c b/src/cell.c
index d4e8f9ff91490cd3a9dd9f78a3e72b19c55bf1cc..2f717850ac423d04ae2d39461c49aa864833d670 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -629,7 +629,8 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, struct cell_buff *buff,
}
/* Re-link the gparts. */
- if (count > 0 && gcount > 0) part_relink_gparts(parts, count, parts_offset);
+ if (count > 0 && gcount > 0)
+ part_relink_gparts_to_parts(parts, count, parts_offset);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that the buffs are OK. */
@@ -744,7 +745,7 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, struct cell_buff *buff,
/* Re-link the parts. */
if (count > 0 && gcount > 0)
- part_relink_parts(gparts, gcount, parts - parts_offset);
+ part_relink_parts_to_gparts(gparts, gcount, parts - parts_offset);
}
/**
diff --git a/src/cell.h b/src/cell.h
index bda0c218975c4b7f534fc50338bfdd2997afc6be..a318cdde2ed95f472c47466dafceff561896588b 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -118,7 +118,7 @@ struct cell {
struct gpart *gparts;
/*! Pointer to the #spart data. */
- struct gpart *sparts;
+ struct spart *sparts;
/*! Pointer for the sorted indices. */
struct entry *sort;
diff --git a/src/common_io.c b/src/common_io.c
index a9981e3024d76c445a8d69b03512a67072296c05..b74d056e98c6bf02bdac80fca1cdab3e9355bde5 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -390,6 +390,8 @@ void writeCodeDescription(hid_t h_file) {
H5Gclose(h_grpcode);
}
+#endif /* HAVE_HDF5 */
+
/* ------------------------------------------------------------------------------------------------
* This part writes the XMF file descriptor enabling a visualisation through
* ParaView
@@ -586,6 +588,9 @@ void prepare_dm_gparts(struct gpart* const gparts, size_t Ndm) {
if (gparts[i].id_or_neg_offset <= 0)
error("0 or negative ID for DM particle %zu: ID=%lld", i,
gparts[i].id_or_neg_offset);
+
+ /* Set gpart type */
+ gparts[i + Ndm].type = swift_type_dark_matter;
}
}
@@ -618,6 +623,9 @@ void duplicate_hydro_gparts(struct part* const parts,
gparts[i + Ndm].mass = hydro_get_mass(&parts[i]);
+ /* Set gpart type */
+ gparts[i + Ndm].type = swift_type_gas;
+
/* Link the particles */
gparts[i + Ndm].id_or_neg_offset = -i;
parts[i].gpart = &gparts[i + Ndm];
@@ -653,6 +661,9 @@ void duplicate_star_gparts(struct spart* const sparts,
gparts[i + Ndm].mass = sparts[i].mass;
+ /* Set gpart type */
+ gparts[i + Ndm].type = swift_type_star;
+
/* Link the particles */
gparts[i + Ndm].id_or_neg_offset = -i;
sparts[i].gpart = &gparts[i + Ndm];
@@ -690,5 +701,3 @@ void collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
error("Collected the wrong number of dm particles (%zu vs. %zu expected)",
count, Ndm);
}
-
-#endif
diff --git a/src/engine.c b/src/engine.c
index 37b6c45598a0f56ccf4ec10fe7717ea357e90cfd..e98ee5090d691fee8a7daadcb1b2247a2b1281d3 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2840,7 +2840,7 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
/* Re-link the gparts. */
if (s->nr_parts > 0 && s->nr_gparts > 0)
- part_relink_gparts(s->parts, s->nr_parts, 0);
+ part_relink_gparts_to_parts(s->parts, s->nr_parts, 0);
/* Re-allocate the local gparts. */
if (e->verbose)
@@ -2857,7 +2857,7 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
/* Re-link the parts. */
if (s->nr_parts > 0 && s->nr_gparts > 0)
- part_relink_parts(s->gparts, s->nr_gparts, s->parts);
+ part_relink_parts_to_gparts(s->gparts, s->nr_gparts, s->parts);
#ifdef SWIFT_DEBUG_CHECKS
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index f06e65e5b30ebcd609c0c6204de33da17b770add..b3953d0e54118ff9fed99cffbc3574eeb122e989 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -25,6 +25,10 @@
/* Gravity particle. */
struct gpart {
+ /* Particle ID. If negative, it is the negative offset of the #part with
+ which this gpart is linked. */
+ long long id_or_neg_offset;
+
/* Particle position. */
double x[3];
@@ -49,9 +53,8 @@ struct gpart {
/* Particle time of end of time-step. */
int ti_end;
- /* Particle ID. If negative, it is the negative offset of the #part with
- which this gpart is linked. */
- long long id_or_neg_offset;
+ /* Type of the #gpart (DM, gas, star, ...) */
+ enum part_type type;
} SWIFT_STRUCT_ALIGN;
diff --git a/src/part.c b/src/part.c
index b00eaccaae0e86f7c4e8019a307f0bf455687b7c..3e0377d3de29868b0c828f0c9ccefdcfbb5b42f5 100644
--- a/src/part.c
+++ b/src/part.c
@@ -36,7 +36,8 @@
* @param N The number of particles to re-link;
* @param offset The offset of #part%s relative to the global parts list.
*/
-void part_relink_gparts(struct part *parts, size_t N, ptrdiff_t offset) {
+void part_relink_gparts_to_parts(struct part *parts, size_t N,
+ ptrdiff_t offset) {
for (size_t k = 0; k < N; k++) {
if (parts[k].gpart) {
parts[k].gpart->id_or_neg_offset = -(k + offset);
@@ -45,20 +46,53 @@ void part_relink_gparts(struct part *parts, size_t N, ptrdiff_t offset) {
}
/**
- * @brief Re-link the #gpart%s associated with the list of #part%s.
+ * @brief Re-link the #gpart%s associated with the list of #spart%s.
+ *
+ * @param sparts The list of #spart.
+ * @param N The number of s-particles to re-link;
+ * @param offset The offset of #spart%s relative to the global sparts list.
+ */
+void part_relink_gparts_to_sparts(struct spart *sparts, size_t N,
+ ptrdiff_t offset) {
+ for (size_t k = 0; k < N; k++) {
+ if (sparts[k].gpart) {
+ sparts[k].gpart->id_or_neg_offset = -(k + offset);
+ }
+ }
+}
+
+/**
+ * @brief Re-link the #part%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 parts The global #part array in which to find the #gpart offsets.
*/
-void part_relink_parts(struct gpart *gparts, size_t N, struct part *parts) {
+void part_relink_parts_to_gparts(struct gpart *gparts, size_t N,
+ struct part *parts) {
for (size_t k = 0; k < N; k++) {
- if (gparts[k].id_or_neg_offset <= 0) {
+ if (gparts[k].type == swift_type_gas) {
parts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
}
}
}
+/**
+ * @brief Re-link the #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 sparts The global #spart array in which to find the #gpart offsets.
+ */
+void part_relink_sparts_to_gparts(struct gpart *gparts, size_t N,
+ struct spart *sparts) {
+ for (size_t k = 0; k < N; k++) {
+ if (gparts[k].type == swift_type_star) {
+ sparts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
+ }
+ }
+}
+
#ifdef WITH_MPI
/* MPI data type for the particle transfers */
MPI_Datatype part_mpi_type;
diff --git a/src/part.h b/src/part.h
index 269e12a9c6b4adbd1046f1a008f1b21ab3214128..149b06eee885cd2614541f5e2afb2db1b8928266 100644
--- a/src/part.h
+++ b/src/part.h
@@ -32,6 +32,7 @@
/* Local headers. */
#include "align.h"
+#include "part_type.h"
/* Some constants. */
#define part_align 128
@@ -66,8 +67,14 @@
/* Import the right star particle definition */
#include "./stars/Default/star_part.h"
-void part_relink_gparts(struct part *parts, size_t N, ptrdiff_t offset);
-void part_relink_parts(struct gpart *gparts, size_t N, struct part *parts);
+void part_relink_gparts_to_parts(struct part *parts, size_t N,
+ ptrdiff_t offset);
+void part_relink_gparts_to_sparts(struct spart *sparts, size_t N,
+ ptrdiff_t offset);
+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);
#ifdef WITH_MPI
/* MPI data type for the particle transfers */
extern MPI_Datatype part_mpi_type;
diff --git a/src/part_type.h b/src/part_type.h
new file mode 100644
index 0000000000000000000000000000000000000000..2c564d6908c8887e8fa8a5197a0a92ed85cbe5bb
--- /dev/null
+++ b/src/part_type.h
@@ -0,0 +1,34 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PART_TYPES_H
+#define SWIFT_PART_TYPES_H
+
+/**
+ * @brief The different types of particles a #gpart can link to.
+ *
+ * Note we use the historical values from Gadget for these fields.
+ */
+enum part_type {
+ swift_type_gas = 0,
+ swift_type_dark_matter = 1,
+ swift_type_star = 4,
+ swift_type_black_hole = 5
+} __attribute__((packed));
+
+#endif /* SWIFT_PART_TYPES_H */
diff --git a/src/space.c b/src/space.c
index 84af59d30fdebe839212ea9119aec92d2384802f..38e72f2e1946e850500c8003f94a508dca47bc9e 100644
--- a/src/space.c
+++ b/src/space.c
@@ -108,6 +108,7 @@ struct parallel_sort {
struct part *parts;
struct gpart *gparts;
struct xpart *xparts;
+ struct spart *sparts;
int *ind;
struct qstack *stack;
unsigned int stack_size;
@@ -211,6 +212,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->sorted = 0;
c->count = 0;
c->gcount = 0;
+ c->scount = 0;
c->init = NULL;
c->extra_ghost = NULL;
c->ghost = NULL;
@@ -389,6 +391,7 @@ void space_regrid(struct space *s, int verbose) {
c->depth = 0;
c->count = 0;
c->gcount = 0;
+ c->scount = 0;
c->super = c;
c->ti_old = ti_current;
lock_init(&c->lock);
@@ -470,6 +473,7 @@ void space_rebuild(struct space *s, int verbose) {
size_t nr_parts = s->nr_parts;
size_t nr_gparts = s->nr_gparts;
+ size_t nr_sparts = s->nr_sparts;
struct cell *restrict cells_top = s->cells_top;
const int ti_current = (s->e != NULL) ? s->e->ti_current : 0;
@@ -490,6 +494,14 @@ void space_rebuild(struct space *s, int verbose) {
if (s->size_gparts > 0)
space_gparts_get_cell_index(s, gind, 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.");
+ if (s->size_sparts > 0)
+ space_sparts_get_cell_index(s, sind, cells_top, verbose);
+
#ifdef WITH_MPI
/* Move non-local parts to the end of the list. */
@@ -612,8 +624,15 @@ void space_rebuild(struct space *s, int verbose) {
if (nr_parts > 0)
space_parts_sort(s, ind, nr_parts, 0, s->nr_cells - 1, verbose);
- /* Re-link the gparts. */
- if (nr_parts > 0 && nr_gparts > 0) part_relink_gparts(s->parts, nr_parts, 0);
+ /* Sort the sparts according to their cells. */
+ if (nr_sparts > 0)
+ space_sparts_sort(s, sind, nr_sparts, 0, s->nr_cells - 1, verbose);
+
+ /* 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);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify space_sort_struct. */
@@ -638,8 +657,19 @@ void space_rebuild(struct space *s, int verbose) {
}
}
+ /* Extract the cell counts from the sorted indices. */
+ size_t last_sindex = 0;
+ sind[nr_parts] = s->nr_cells; // sentinel.
+ for (size_t k = 0; k < nr_sparts; k++) {
+ if (sind[k] < sind[k + 1]) {
+ cells_top[sind[k]].scount = k - last_sindex + 1;
+ last_sindex = k + 1;
+ }
+ }
+
/* We no longer need the indices as of here. */
free(ind);
+ free(sind);
#ifdef WITH_MPI
@@ -667,7 +697,7 @@ void space_rebuild(struct space *s, int verbose) {
}
nr_gparts = s->nr_gparts;
-#endif
+#endif /* WITH_MPI */
/* Sort the gparts according to their cells. */
if (nr_gparts > 0)
@@ -675,7 +705,11 @@ void space_rebuild(struct space *s, int verbose) {
/* Re-link the parts. */
if (nr_parts > 0 && nr_gparts > 0)
- part_relink_parts(s->gparts, nr_gparts, s->parts);
+ 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;
@@ -719,15 +753,18 @@ void space_rebuild(struct space *s, int verbose) {
struct part *finger = s->parts;
struct xpart *xfinger = s->xparts;
struct gpart *gfinger = s->gparts;
+ struct spart *sfinger = s->sparts;
for (int k = 0; k < s->nr_cells; k++) {
struct cell *restrict c = &cells_top[k];
c->ti_old = ti_current;
c->parts = finger;
c->xparts = xfinger;
c->gparts = gfinger;
+ c->sparts = sfinger;
finger = &finger[c->count];
xfinger = &xfinger[c->count];
gfinger = &gfinger[c->gcount];
+ sfinger = &sfinger[c->scount];
}
// message( "hooking up cells took %.3f %s." ,
// clocks_from_ticks(getticks() - tic), clocks_getunit());
@@ -892,8 +929,58 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts,
}
/**
- * @brief Computes the cell index of all the particles and update the cell
- * count.
+ * @brief #threadpool mapper function to compute the s-particle cell indices.
+ *
+ * @param map_data Pointer towards the s-particles.
+ * @param nr_sparts The number of s-particles to treat.
+ * @param extra_data Pointers to the space and index list
+ */
+void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
+ void *extra_data) {
+
+ /* Unpack the data */
+ struct spart *restrict sparts = (struct spart *)map_data;
+ struct index_data *data = (struct index_data *)extra_data;
+ struct space *s = data->s;
+ int *const ind = data->ind + (ptrdiff_t)(sparts - s->sparts);
+
+ /* Get some constants */
+ const double dim_x = s->dim[0];
+ const double dim_y = s->dim[1];
+ const double dim_z = s->dim[2];
+ const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
+ const double ih_x = s->iwidth[0];
+ const double ih_y = s->iwidth[1];
+ const double ih_z = s->iwidth[2];
+
+ for (int k = 0; k < nr_sparts; k++) {
+
+ /* Get the particle */
+ struct spart *restrict sp = &sparts[k];
+
+ const double old_pos_x = sp->x[0];
+ const double old_pos_y = sp->x[1];
+ const double old_pos_z = sp->x[2];
+
+ /* Put it back into the simulation volume */
+ const double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
+ const double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
+ const double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+
+ /* Get its cell index */
+ const int index =
+ cell_getid(cdim, pos_x * ih_x, pos_y * ih_y, pos_z * ih_z);
+ ind[k] = index;
+
+ /* Update the position */
+ sp->x[0] = pos_x;
+ sp->x[1] = pos_y;
+ sp->x[2] = pos_z;
+ }
+}
+
+/**
+ * @brief Computes the cell index of all the particles.
*
* @param s The #space.
* @param ind The array of indices to fill.
@@ -920,8 +1007,7 @@ void space_parts_get_cell_index(struct space *s, int *ind, struct cell *cells,
}
/**
- * @brief Computes the cell index of all the g-particles and update the cell
- * gcount.
+ * @brief Computes the cell index of all the g-particles.
*
* @param s The #space.
* @param gind The array of indices to fill.
@@ -947,6 +1033,33 @@ void space_gparts_get_cell_index(struct space *s, int *gind, struct cell *cells,
clocks_getunit());
}
+/**
+ * @brief Computes the cell index of all the s-particles.
+ *
+ * @param s The #space.
+ * @param sind The array of indices to fill.
+ * @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) {
+
+ const ticks tic = getticks();
+
+ /* Pack the extra information */
+ struct index_data data;
+ data.s = s;
+ data.cells = cells;
+ data.ind = sind;
+
+ threadpool_map(&s->e->threadpool, space_sparts_get_cell_index_mapper,
+ s->sparts, s->nr_sparts, sizeof(struct gpart), 1000, &data);
+
+ if (verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+}
+
/**
* @brief Sort the particles and condensed particles according to the given
* indices.
@@ -1127,6 +1240,182 @@ void space_parts_sort_mapper(void *map_data, int num_elements,
} /* main loop. */
}
+/**
+ * @brief Sort the s-particles according to the given indices.
+ *
+ * @param s The #space.
+ * @param ind The indices with respect to which the #spart are sorted.
+ * @param N The number of parts
+ * @param min Lowest index.
+ * @param max highest index.
+ * @param verbose Are we talkative ?
+ */
+void space_sparts_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.sparts = s->sparts;
+ sort_struct.ind = ind;
+ sort_struct.stack_size = 2 * (max - min + 1) + 10 + s->e->nr_threads;
+ if ((sort_struct.stack =
+ 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_sparts_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);
+ 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_sparts_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 spart *sparts = sort_struct->sparts;
+
+ /* 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) {
+ memswap(&ind[ii], &ind[jj], sizeof(int));
+ memswap(&sparts[ii], &sparts[jj], sizeof(struct spart));
+ }
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Verify space_sort_struct. */
+ for (int k = i; k <= jj; k++)
+ if (ind[k] > pivot) {
+ message("sorting failed at k=%i, ind[k]=%i, pivot=%i, i=%li, j=%li.",
+ k, ind[k], pivot, i, j);
+ error("Partition failed (<=pivot).");
+ }
+ for (int k = jj + 1; k <= j; k++)
+ if (ind[k] <= pivot) {
+ message("sorting failed at k=%i, ind[k]=%i, pivot=%i, i=%li, j=%li.",
+ k, ind[k], pivot, i, j);
+ error("Partition failed (>pivot).");
+ }
+#endif
+
+ /* Split-off largest interval. */
+ if (jj - i > j - jj + 1) {
+
+ /* Recurse on the left? */
+ if (jj > i && pivot > min) {
+ qid = atomic_inc(&sort_struct->last) % sort_struct->stack_size;
+ while (sort_struct->stack[qid].ready)
+ ;
+ sort_struct->stack[qid].i = i;
+ sort_struct->stack[qid].j = jj;
+ sort_struct->stack[qid].min = min;
+ sort_struct->stack[qid].max = pivot;
+ if (atomic_inc(&sort_struct->waiting) >= sort_struct->stack_size)
+ error("Qstack overflow.");
+ sort_struct->stack[qid].ready = 1;
+ }
+
+ /* Recurse on the right? */
+ if (jj + 1 < j && pivot + 1 < max) {
+ i = jj + 1;
+ min = pivot + 1;
+ } else
+ break;
+
+ } else {
+
+ /* Recurse on the right? */
+ if (pivot + 1 < max) {
+ qid = atomic_inc(&sort_struct->last) % sort_struct->stack_size;
+ while (sort_struct->stack[qid].ready)
+ ;
+ sort_struct->stack[qid].i = jj + 1;
+ sort_struct->stack[qid].j = j;
+ sort_struct->stack[qid].min = pivot + 1;
+ sort_struct->stack[qid].max = max;
+ if (atomic_inc(&sort_struct->waiting) >= sort_struct->stack_size)
+ error("Qstack overflow.");
+ sort_struct->stack[qid].ready = 1;
+ }
+
+ /* Recurse on the left? */
+ if (jj > i && pivot > min) {
+ j = jj;
+ max = pivot;
+ } else
+ break;
+ }
+
+ } /* loop over sub-intervals. */
+
+ atomic_dec(&sort_struct->waiting);
+
+ } /* main loop. */
+}
+
/**
* @brief Sort the g-particles according to the given indices.
*
diff --git a/src/space.h b/src/space.h
index 1103cc4867c0ff7ded700331a1a4ff540233e8a7..17e4fc507488ba9f494dbf8d6ca0b702b2bc543e 100644
--- a/src/space.h
+++ b/src/space.h
@@ -153,6 +153,8 @@ void space_parts_sort(struct space *s, int *ind, size_t N, int min, int max,
int verbose);
void space_gparts_sort(struct space *s, int *ind, size_t N, int min, int max,
int verbose);
+void space_sparts_sort(struct space *s, int *ind, size_t N, int min, int max,
+ int verbose);
struct cell *space_getcell(struct space *s);
int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
double *shift);
@@ -176,6 +178,8 @@ void space_parts_sort_mapper(void *map_data, int num_elements,
void *extra_data);
void space_gparts_sort_mapper(void *map_data, int num_elements,
void *extra_data);
+void space_sparts_sort_mapper(void *map_data, int num_elements,
+ void *extra_data);
void space_rebuild(struct space *s, int verbose);
void space_recycle(struct space *s, struct cell *c);
void space_recycle_list(struct space *s, struct cell *list_begin,
@@ -187,8 +191,11 @@ void space_parts_get_cell_index(struct space *s, int *ind, struct cell *cells,
int verbose);
void space_gparts_get_cell_index(struct space *s, int *gind, struct cell *cells,
int verbose);
+void space_sparts_get_cell_index(struct space *s, int *sind, struct cell *cells,
+ int verbose);
void space_do_parts_sort();
void space_do_gparts_sort();
+void space_do_sparts_sort();
void space_init_parts(struct space *s);
void space_init_gparts(struct space *s);
void space_init_sparts(struct space *s);
diff --git a/src/stars/Default/star_part.h b/src/stars/Default/star_part.h
index 3ee5a9452549f5f02aaf8131204b5018ccb1d029..943a392db2e6e93a52c95f4aa4ec794246589939 100644
--- a/src/stars/Default/star_part.h
+++ b/src/stars/Default/star_part.h
@@ -40,6 +40,7 @@ struct spart {
/* Particle velocity. */
float v[3];
+ /* Star mass */
float mass;
/* Particle time of beginning of time-step. */