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

Star particles are now sorted among top-level cells and linked correctly to their gpart

parent 47e86b44
......@@ -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 \
......
......@@ -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);
}
/**
......
......@@ -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;
......
......@@ -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
......@@ -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
......
......@@ -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;
......
......@@ -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;
......
......@@ -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;
......
/*******************************************************************************
* 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 */
......@@ -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 {