Commit 8a86e252 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' into parallel_gpart_sort

parents c2790396 31fe964e
......@@ -363,7 +363,7 @@ int main(int argc, char *argv[]) {
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
read_ic_parallel(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
#else
read_ic_serial(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
......@@ -396,8 +396,8 @@ int main(int argc, char *argv[]) {
/* MATTHIEU: Temporary fix to preserve master */
if (!with_gravity) {
free(gparts);
for(size_t k = 0; k < Ngas; ++k)
parts[k].gpart = NULL;
gparts = NULL;
for (size_t k = 0; k < Ngas; ++k) parts[k].gpart = NULL;
Ngpart = 0;
#if defined(WITH_MPI)
N_long[0] = Ngas;
......
......@@ -89,14 +89,18 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
c->ti_end_min = pc->ti_end_min;
c->ti_end_max = pc->ti_end_max;
c->count = pc->count;
c->gcount = pc->gcount;
c->tag = pc->tag;
/* Fill the progeny recursively, depth-first. */
/* Number of new cells created. */
int count = 1;
/* Fill the progeny recursively, depth-first. */
for (int k = 0; k < 8; k++)
if (pc->progeny[k] >= 0) {
struct cell *temp = space_getcell(s);
temp->count = 0;
temp->gcount = 0;
temp->loc[0] = c->loc[0];
temp->loc[1] = c->loc[1];
temp->loc[2] = c->loc[2];
......@@ -122,7 +126,7 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
}
/**
* @brief Link the cells recursively to the given part array.
* @brief Link the cells recursively to the given #part array.
*
* @param c The #cell.
* @param parts The #part array.
......@@ -130,7 +134,7 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
* @return The number of particles linked.
*/
int cell_link(struct cell *c, struct part *parts) {
int cell_link_parts(struct cell *c, struct part *parts) {
c->parts = parts;
......@@ -139,14 +143,40 @@ int cell_link(struct cell *c, struct part *parts) {
int offset = 0;
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL)
offset += cell_link(c->progeny[k], &parts[offset]);
offset += cell_link_parts(c->progeny[k], &parts[offset]);
}
}
/* Return the total number of unpacked cells. */
/* Return the total number of linked particles. */
return c->count;
}
/**
* @brief Link the cells recursively to the given #gpart array.
*
* @param c The #cell.
* @param gparts The #gpart array.
*
* @return The number of particles linked.
*/
int cell_link_gparts(struct cell *c, struct gpart *gparts) {
c->gparts = gparts;
/* Fill the progeny recursively, depth-first. */
if (c->split) {
int offset = 0;
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL)
offset += cell_link_gparts(c->progeny[k], &gparts[offset]);
}
}
/* Return the total number of linked particles. */
return c->gcount;
}
/**
* @brief Pack the data of the given cell and all it's sub-cells.
*
......@@ -164,6 +194,7 @@ int cell_pack(struct cell *c, struct pcell *pc) {
pc->ti_end_min = c->ti_end_min;
pc->ti_end_max = c->ti_end_max;
pc->count = c->count;
pc->gcount = c->gcount;
c->tag = pc->tag = atomic_inc(&cell_next_tag) % cell_max_tag;
/* Fill in the progeny, depth-first recursion. */
......
......@@ -44,7 +44,7 @@ struct pcell {
int ti_end_min, ti_end_max;
/* Number of particles in this cell. */
int count;
int count, gcount;
/* tag used for MPI communication. */
int tag;
......@@ -175,7 +175,8 @@ void cell_gunlocktree(struct cell *c);
int cell_pack(struct cell *c, struct pcell *pc);
int cell_unpack(struct pcell *pc, struct cell *c, struct space *s);
int cell_getsize(struct cell *c);
int cell_link(struct cell *c, struct part *parts);
int cell_link_parts(struct cell *c, struct part *parts);
int cell_link_gparts(struct cell *c, struct gpart *gparts);
void cell_init_parts(struct cell *c, void *data);
void cell_convert_hydro(struct cell *c, void *data);
void cell_clean_links(struct cell *c, void *data);
......
......@@ -502,7 +502,7 @@ void writeXMFline(FILE* xmfFile, char* fileName, char* partTypeGroupName,
* @param gparts The array of #gpart freshly read in.
* @param Ndm The number of DM particles read in.
*/
void prepare_dm_gparts(struct gpart* gparts, size_t Ndm) {
void prepare_dm_gparts(struct gpart* const gparts, size_t Ndm) {
/* Let's give all these gparts a negative id */
for (size_t i = 0; i < Ndm; ++i) {
......@@ -527,8 +527,9 @@ void prepare_dm_gparts(struct gpart* gparts, size_t Ndm) {
* @param Ngas The number of gas particles read in.
* @param Ndm The number of DM particles read in.
*/
void duplicate_hydro_gparts(struct part* parts, struct gpart* gparts,
size_t Ngas, size_t Ndm) {
void duplicate_hydro_gparts(struct part* const parts,
struct gpart* const gparts, size_t Ngas,
size_t Ndm) {
for (size_t i = 0; i < Ngas; ++i) {
......@@ -557,16 +558,19 @@ void duplicate_hydro_gparts(struct part* parts, struct gpart* gparts,
* @param dmparts The array of #gpart containg DM particles to be filled.
* @param Ndm The number of DM particles.
*/
void collect_dm_gparts(struct gpart* gparts, size_t Ntot, struct gpart* dmparts,
size_t Ndm) {
void collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
struct gpart* const dmparts, size_t Ndm) {
size_t count = 0;
/* Loop over all gparts */
for (size_t i = 0; i < Ntot; ++i) {
/* message("i=%zd count=%zd id=%lld part=%p", i, count, gparts[i].id,
* gparts[i].part); */
/* And collect the DM ones */
if (gparts[i].id < 0) {
if (gparts[i].id < 0LL) {
memcpy(&dmparts[count], &gparts[i], sizeof(struct gpart));
dmparts[count].id = -dmparts[count].id;
count++;
......
......@@ -78,11 +78,12 @@ extern const char* particle_type_names[];
hid_t hdf5Type(enum DATA_TYPE type);
size_t sizeOfType(enum DATA_TYPE type);
void collect_dm_gparts(struct gpart* gparts, size_t Ntot, struct gpart* dmparts,
size_t Ndm);
void prepare_dm_gparts(struct gpart* gparts, size_t Ndm);
void duplicate_hydro_gparts(struct part* parts, struct gpart* gparts,
size_t Ngas, size_t Ndm);
void collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
struct gpart* const dmparts, size_t Ndm);
void prepare_dm_gparts(struct gpart* const gparts, size_t Ndm);
void duplicate_hydro_gparts(struct part* const parts,
struct gpart* const gparts, size_t Ngas,
size_t Ndm);
void readAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data);
......
This diff is collapsed.
......@@ -62,6 +62,8 @@ extern const char *engine_policy_names[];
#define engine_maxtaskspercell 96
#define engine_maxproxies 64
#define engine_tasksreweight 10
#define engine_parts_size_grow 1.05
#define engine_redistribute_alloc_margin 1.2
/* The rank of the engine as a global variable (for messages). */
extern int engine_rank;
......@@ -160,12 +162,6 @@ struct engine {
/* Are we talkative ? */
int verbose;
#ifdef WITH_MPI
/* MPI data type for the particle transfers */
MPI_Datatype part_mpi_type;
MPI_Datatype xpart_mpi_type;
#endif
};
/* Function prototypes. */
......@@ -182,7 +178,9 @@ void engine_init_particles(struct engine *e);
void engine_step(struct engine *e);
void engine_maketasks(struct engine *e);
void engine_split(struct engine *e, struct partition *initial_partition);
int engine_exchange_strays(struct engine *e, int offset, size_t *ind, size_t N);
void engine_exchange_strays(struct engine *e, size_t offset_parts,
int *ind_part, size_t *Npart, size_t offset_gparts,
int *ind_gpart, size_t *Ngpart);
void engine_rebuild(struct engine *e);
void engine_repartition(struct engine *e);
void engine_makeproxies(struct engine *e);
......
......@@ -26,33 +26,21 @@
#endif
/* This object's header. */
#include "error.h"
#include "part.h"
#ifdef WITH_MPI
/**
* @brief Registers and returns an MPI type for the particles
*
* @param part_type The type container
*/
void part_create_mpi_type(MPI_Datatype* part_type) {
/* This is not the recommended way of doing this.
One should define the structure field by field
But as long as we don't do serialization via MPI-IO
we don't really care.
Also we would have to modify this function everytime something
is added to the part structure. */
MPI_Type_contiguous(sizeof(struct part) / sizeof(unsigned char), MPI_BYTE,
part_type);
MPI_Type_commit(part_type);
}
/* MPI data type for the particle transfers */
MPI_Datatype part_mpi_type;
MPI_Datatype xpart_mpi_type;
MPI_Datatype gpart_mpi_type;
#endif
#ifdef WITH_MPI
/**
* @brief Registers and returns an MPI type for the xparticles
*
* @param xpart_type The type container
* @brief Registers MPI particle types.
*/
void xpart_create_mpi_type(MPI_Datatype* xpart_type) {
void part_create_mpi_types() {
/* This is not the recommended way of doing this.
One should define the structure field by field
......@@ -60,9 +48,20 @@ void xpart_create_mpi_type(MPI_Datatype* xpart_type) {
we don't really care.
Also we would have to modify this function everytime something
is added to the part structure. */
MPI_Type_contiguous(sizeof(struct xpart) / sizeof(unsigned char), MPI_BYTE,
xpart_type);
MPI_Type_commit(xpart_type);
if (MPI_Type_contiguous(sizeof(struct part) / sizeof(unsigned char), MPI_BYTE,
&part_mpi_type) != MPI_SUCCESS ||
MPI_Type_commit(&part_mpi_type) != MPI_SUCCESS) {
error("Failed to create MPI type for parts.");
}
if (MPI_Type_contiguous(sizeof(struct xpart) / sizeof(unsigned char),
MPI_BYTE, &xpart_mpi_type) != MPI_SUCCESS ||
MPI_Type_commit(&xpart_mpi_type) != MPI_SUCCESS) {
error("Failed to create MPI type for xparts.");
}
if (MPI_Type_contiguous(sizeof(struct gpart) / sizeof(unsigned char),
MPI_BYTE, &gpart_mpi_type) != MPI_SUCCESS ||
MPI_Type_commit(&gpart_mpi_type) != MPI_SUCCESS) {
error("Failed to create MPI type for gparts.");
}
}
#endif
......@@ -35,8 +35,8 @@
/* Some constants. */
#define part_align 64
#define gpart_align 32
#define xpart_align 32
#define gpart_align 32
/* Import the right particle definition */
#if defined(MINIMAL_SPH)
......@@ -52,8 +52,12 @@
#include "./gravity/Default/gravity_part.h"
#ifdef WITH_MPI
void part_create_mpi_type(MPI_Datatype* part_type);
void xpart_create_mpi_type(MPI_Datatype* xpart_type);
/* MPI data type for the particle transfers */
extern MPI_Datatype part_mpi_type;
extern MPI_Datatype xpart_mpi_type;
extern MPI_Datatype gpart_mpi_type;
void part_create_mpi_types();
#endif
#endif /* SWIFT_PART_H */
......@@ -189,20 +189,21 @@ void proxy_parts_exch1(struct proxy *p) {
#ifdef WITH_MPI
/* Send the number of particles. */
if (MPI_Isend(&p->nr_parts_out, 1, MPI_INT, p->nodeID,
p->buff_out[0] = p->nr_parts_out;
p->buff_out[1] = p->nr_gparts_out;
if (MPI_Isend(p->buff_out, 2, MPI_INT, p->nodeID,
p->mynodeID * proxy_tag_shift + proxy_tag_count, MPI_COMM_WORLD,
&p->req_parts_count_out) != MPI_SUCCESS)
error("Failed to isend nr of parts.");
// message( "isent particle count (%i) from node %i to node %i." ,
// p->nr_parts_out , p->mynodeID , p->nodeID ); fflush(stdout);
/* message( "isent particle counts [%i, %i] from node %i to node %i." ,
p->buff_out[0], p->buff_out[1], p->mynodeID , p->nodeID ); fflush(stdout); */
/* Send the particle buffers. */
if (p->nr_parts_out > 0) {
if (MPI_Isend(p->parts_out, sizeof(struct part) * p->nr_parts_out, MPI_BYTE,
p->nodeID, p->mynodeID * proxy_tag_shift + proxy_tag_parts,
if (MPI_Isend(p->parts_out, p->nr_parts_out, part_mpi_type, p->nodeID,
p->mynodeID * proxy_tag_shift + proxy_tag_parts,
MPI_COMM_WORLD, &p->req_parts_out) != MPI_SUCCESS ||
MPI_Isend(p->xparts_out, sizeof(struct xpart) * p->nr_parts_out,
MPI_BYTE, p->nodeID,
MPI_Isend(p->xparts_out, p->nr_parts_out, xpart_mpi_type, p->nodeID,
p->mynodeID * proxy_tag_shift + proxy_tag_xparts,
MPI_COMM_WORLD, &p->req_xparts_out) != MPI_SUCCESS)
error("Failed to isend part data.");
......@@ -213,14 +214,20 @@ void proxy_parts_exch1(struct proxy *p) {
p->parts_out[k].id, p->parts_out[k].x[0], p->parts_out[k].x[1],
p->parts_out[k].x[2], p->parts_out[k].h, p->nodeID);*/
}
if (p->nr_gparts_out > 0) {
if (MPI_Isend(p->gparts_out, p->nr_gparts_out, gpart_mpi_type, p->nodeID,
p->mynodeID * proxy_tag_shift + proxy_tag_gparts,
MPI_COMM_WORLD, &p->req_gparts_out) != MPI_SUCCESS)
error("Failed to isend part data.");
// message( "isent gpart data (%i) to node %i." , p->nr_parts_out ,
// p->nodeID ); fflush(stdout);
}
/* Receive the number of particles. */
if (MPI_Irecv(&p->nr_parts_in, 1, MPI_INT, p->nodeID,
if (MPI_Irecv(p->buff_in, 2, MPI_INT, p->nodeID,
p->nodeID * proxy_tag_shift + proxy_tag_count, MPI_COMM_WORLD,
&p->req_parts_count_in) != MPI_SUCCESS)
error("Failed to irecv nr of parts.");
// message( "irecv particle count on node %i from node %i." , p->mynodeID ,
// p->nodeID ); fflush(stdout);
#else
error("SWIFT was not compiled with MPI support.");
......@@ -231,6 +238,10 @@ void proxy_parts_exch2(struct proxy *p) {
#ifdef WITH_MPI
/* Unpack the incomming parts counts. */
p->nr_parts_in = p->buff_in[0];
p->nr_gparts_in = p->buff_in[1];
/* Is there enough space in the buffer? */
if (p->nr_parts_in > p->size_parts_in) {
do {
......@@ -244,19 +255,36 @@ void proxy_parts_exch2(struct proxy *p) {
p->size_parts_in)) == NULL)
error("Failed to re-allocate parts_in buffers.");
}
if (p->nr_gparts_in > p->size_gparts_in) {
do {
p->size_gparts_in *= proxy_buffgrow;
} while (p->nr_gparts_in > p->size_gparts_in);
free(p->gparts_in);
if ((p->gparts_in = (struct gpart *)malloc(sizeof(struct gpart) *
p->size_gparts_in)) == NULL)
error("Failed to re-allocate gparts_in buffers.");
}
/* Receive the particle buffers. */
if (p->nr_parts_in > 0) {
if (MPI_Irecv(p->parts_in, sizeof(struct part) * p->nr_parts_in, MPI_BYTE,
p->nodeID, p->nodeID * proxy_tag_shift + proxy_tag_parts,
MPI_COMM_WORLD, &p->req_parts_in) != MPI_SUCCESS ||
MPI_Irecv(p->xparts_in, sizeof(struct xpart) * p->nr_parts_in, MPI_BYTE,
p->nodeID, p->nodeID * proxy_tag_shift + proxy_tag_xparts,
if (MPI_Irecv(p->parts_in, p->nr_parts_in, part_mpi_type, p->nodeID,
p->nodeID * proxy_tag_shift + proxy_tag_parts, MPI_COMM_WORLD,
&p->req_parts_in) != MPI_SUCCESS ||
MPI_Irecv(p->xparts_in, p->nr_parts_in, xpart_mpi_type, p->nodeID,
p->nodeID * proxy_tag_shift + proxy_tag_xparts,
MPI_COMM_WORLD, &p->req_xparts_in) != MPI_SUCCESS)
error("Failed to irecv part data.");
// message( "irecv particle data (%i) from node %i." , p->nr_parts_in ,
// p->nodeID ); fflush(stdout);
}
if (p->nr_gparts_in > 0) {
if (MPI_Irecv(p->gparts_in, p->nr_gparts_in, gpart_mpi_type, p->nodeID,
p->nodeID * proxy_tag_shift + proxy_tag_gparts,
MPI_COMM_WORLD, &p->req_gparts_in) != MPI_SUCCESS)
error("Failed to irecv gpart data.");
// message( "irecv gpart data (%i) from node %i." , p->nr_gparts_in ,
// p->nodeID ); fflush(stdout);
}
#else
error("SWIFT was not compiled with MPI support.");
......@@ -303,6 +331,37 @@ void proxy_parts_load(struct proxy *p, const struct part *parts,
p->nr_parts_out += N;
}
/**
* @brief Load parts onto a proxy for exchange.
*
* @param p The #proxy.
* @param gparts Pointer to an array of #gpart to send.
* @param N The number of parts.
*/
void proxy_gparts_load(struct proxy *p, const struct gpart *gparts, int N) {
/* Is there enough space in the buffer? */
if (p->nr_gparts_out + N > p->size_gparts_out) {
do {
p->size_gparts_out *= proxy_buffgrow;
} while (p->nr_gparts_out + N > p->size_gparts_out);
struct gpart *tp;
if ((tp = (struct gpart *)malloc(sizeof(struct gpart) *
p->size_gparts_out)) == NULL)
error("Failed to re-allocate gparts_out buffers.");
memcpy(tp, p->gparts_out, sizeof(struct gpart) * p->nr_gparts_out);
free(p->gparts_out);
p->gparts_out = tp;
}
/* Copy the parts and xparts data to the buffer. */
memcpy(&p->gparts_out[p->nr_gparts_out], gparts, sizeof(struct gpart) * N);
/* Increase the counters. */
p->nr_gparts_out += N;
}
/**
* @brief Initialize the given proxy.
*
......@@ -352,4 +411,20 @@ void proxy_init(struct proxy *p, int mynodeID, int nodeID) {
error("Failed to allocate parts_out buffers.");
}
p->nr_parts_out = 0;
/* Allocate the gpart send and receive buffers, if needed. */
if (p->gparts_in == NULL) {
p->size_gparts_in = proxy_buffinit;
if ((p->gparts_in = (struct gpart *)malloc(sizeof(struct gpart) *
p->size_gparts_in)) == NULL)
error("Failed to allocate gparts_in buffers.");
}
p->nr_gparts_in = 0;
if (p->gparts_out == NULL) {
p->size_gparts_out = proxy_buffinit;
if ((p->gparts_out = (struct gpart *)malloc(sizeof(struct gpart) *
p->size_gparts_out)) == NULL)
error("Failed to allocate gparts_out buffers.");
}
p->nr_gparts_out = 0;
}
......@@ -32,7 +32,8 @@
#define proxy_tag_count 0
#define proxy_tag_parts 1
#define proxy_tag_xparts 2
#define proxy_tag_cells 3
#define proxy_tag_gparts 3
#define proxy_tag_cells 4
/* Data structure for the proxy. */
struct proxy {
......@@ -53,14 +54,21 @@ struct proxy {
/* The parts and xparts buffers for input and output. */
struct part *parts_in, *parts_out;
struct xpart *xparts_in, *xparts_out;
struct gpart *gparts_in, *gparts_out;
int size_parts_in, size_parts_out;
int nr_parts_in, nr_parts_out;
int size_gparts_in, size_gparts_out;
int nr_gparts_in, nr_gparts_out;
/* Buffer to hold the incomming/outgoing particle counts. */
int buff_out[2], buff_in[2];
/* MPI request handles. */
#ifdef WITH_MPI
MPI_Request req_parts_count_out, req_parts_count_in;
MPI_Request req_parts_out, req_parts_in;
MPI_Request req_xparts_out, req_xparts_in;
MPI_Request req_gparts_out, req_gparts_in;
MPI_Request req_cells_count_out, req_cells_count_in;
MPI_Request req_cells_out, req_cells_in;
#endif
......@@ -70,6 +78,7 @@ struct proxy {
void proxy_init(struct proxy *p, int mynodeID, int nodeID);
void proxy_parts_load(struct proxy *p, const struct part *parts,
const struct xpart *xparts, int N);
void proxy_gparts_load(struct proxy *p, const struct gpart *gparts, int N);
void proxy_parts_exch1(struct proxy *p);
void proxy_parts_exch2(struct proxy *p);
void proxy_addcell_in(struct proxy *p, struct cell *c);
......
......@@ -1098,7 +1098,7 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
break;
case task_type_recv:
#ifdef WITH_MPI
err = MPI_Irecv(t->ci->parts, t->ci->count, s->part_mpi_type,
err = MPI_Irecv(t->ci->parts, t->ci->count, part_mpi_type,
t->ci->nodeID, t->flags, MPI_COMM_WORLD, &t->req);
if (err != MPI_SUCCESS) {
mpi_error(err, "Failed to emit irecv for particle data.");
......@@ -1113,7 +1113,7 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
break;
case task_type_send:
#ifdef WITH_MPI
err = MPI_Isend(t->ci->parts, t->ci->count, s->part_mpi_type,
err = MPI_Isend(t->ci->parts, t->ci->count, part_mpi_type,
t->cj->nodeID, t->flags, MPI_COMM_WORLD, &t->req);
if (err != MPI_SUCCESS) {
mpi_error(err, "Failed to emit isend for particle data.");
......@@ -1354,12 +1354,6 @@ void scheduler_init(struct scheduler *s, struct space *space, int nr_tasks,
s->tasks = NULL;
s->tasks_ind = NULL;
scheduler_reset(s, nr_tasks);
/* Construct types for MPI communications */
#ifdef WITH_MPI
part_create_mpi_type(&s->part_mpi_type);
xpart_create_mpi_type(&s->xpart_mpi_type);
#endif
}
/**
......
......@@ -100,12 +100,6 @@ struct scheduler {
/* The node we are working on. */
int nodeID;
#ifdef WITH_MPI
/* MPI data type for the particle transfers */
MPI_Datatype part_mpi_type;
MPI_Datatype xpart_mpi_type;
#endif
};
/* Function prototypes. */
......
......@@ -305,7 +305,7 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
*/
void space_rebuild(struct space *s, double cell_max, int verbose) {
const ticks tic = getticks();
/* Be verbose about this. */
......@@ -318,23 +318,15 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
int nr_gparts = s->nr_gparts;
struct cell *restrict cells = s->cells;
double ih[3], dim[3];
int cdim[3];
ih[0] = s->ih[0];
ih[1] = s->ih[1];
ih[2] = s->ih[2];
dim[0] = s->dim[0];
dim[1] = s->dim[1];
dim[2] = s->dim[2];
cdim[0] = s->cdim[0];
cdim[1] = s->cdim[1];
cdim[2] = s->cdim[2];
const double ih[3] = {s->ih[0], s->ih[1], s->ih[2]};
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
/* Run through the particles and get their cell index. */
// tic = getticks();
const size_t ind_size = s->size_parts;
size_t *ind;
if ((ind = (size_t *)malloc(sizeof(size_t) * ind_size)) == NULL)
int *ind;
if ((ind = (int *)malloc(sizeof(int) * ind_size)) == NULL)
error("Failed to allocate temporary particle indices.");
for (int k = 0; k < nr_parts; k++) {
struct part *restrict p = &s->parts[k];
......@@ -347,37 +339,91 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
cells[ind[k]].count++;
}
// message( "getting particle indices took %.3f %s." ,
// clocks_from_ticks(getticks() - tic), clocks_getunit()):
/* Run through the gravity particles and get their cell index. */
// tic = getticks();
const size_t gind_size = s->size_gparts;
int *gind;
if ((gind = (int *)malloc(sizeof(int) * gind_size)) == NULL)
error("Failed to allocate temporary g-particle indices.");
for (int k = 0; k < nr_gparts; k++) {