Commit 7236e2db authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'gpart_ids' into 'master'

One major change to replace all gpart id/parts with the new id_or_neg_offset.

As discussed in #123, this merge request changes the `gpart.id` and `gpart.part` anonymous union into a single `long long int id_or_neg_offset`.

If this number is positive, then it is the particle ID and the `gpart` has no hydro counterpart.

If this number is negative, then it is the offset in the `parts` array of this particle's hydro counterpart.

See merge request !118
parents 6d9f54ae 19a99afb
......@@ -354,7 +354,7 @@ int main(int argc, char *argv[]) {
free(parts);
parts = NULL;
for (size_t k = 0; k < Ngpart; ++k)
if (gparts[k].id > 0) error("Linking problem");
if (gparts[k].id_or_neg_offset < 0) error("Linking problem");
Ngas = 0;
}
......
......@@ -376,9 +376,11 @@ void cell_gunlocktree(struct cell *c) {
* @brief Sort the parts into eight bins along the given pivots.
*
* @param c The #cell array to be sorted.
* @param parts_offset Offset of the cell parts array relative to the
* space's parts array, i.e. c->parts - s->parts.
*/
void cell_split(struct cell *c) {
void cell_split(struct cell *c, ptrdiff_t parts_offset) {
int i, j;
const int count = c->count, gcount = c->gcount;
......@@ -496,8 +498,7 @@ void cell_split(struct cell *c) {
}
/* Re-link the gparts. */
for (int k = 0; k < count; k++)
if (parts[k].gpart != NULL) parts[k].gpart->part = &parts[k];
part_relink_gparts(parts, count, parts_offset);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that _all_ the parts have been assigned to a cell. */
......@@ -592,8 +593,7 @@ void cell_split(struct cell *c) {
}
/* Re-link the parts. */
for (int k = 0; k < gcount; k++)
if (gparts[k].id > 0) gparts[k].part->gpart = &gparts[k];
part_relink_parts(gparts, gcount, parts - parts_offset);
}
/**
......
......@@ -24,6 +24,9 @@
#define SWIFT_CELL_H
/* Includes. */
#include <stddef.h>
/* Local includes. */
#include "lock.h"
#include "multipole.h"
#include "part.h"
......@@ -175,7 +178,7 @@ struct cell {
((int)(k) + (cdim)[2] * ((int)(j) + (cdim)[1] * (int)(i)))
/* Function prototypes. */
void cell_split(struct cell *c);
void cell_split(struct cell *c, ptrdiff_t parts_offset);
int cell_locktree(struct cell *c);
void cell_unlocktree(struct cell *c);
int cell_glocktree(struct cell *c);
......
......@@ -516,12 +516,10 @@ 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) {
/* 0 or negative ids are not allowed */
if (gparts[i].id <= 0)
error("0 or negative ID for DM particle %zd: ID=%lld", i, gparts[i].id);
gparts[i].id = -gparts[i].id;
if (gparts[i].id_or_neg_offset <= 0)
error("0 or negative ID for DM particle %zd: ID=%lld", i,
gparts[i].id_or_neg_offset);
}
}
......@@ -555,7 +553,7 @@ void duplicate_hydro_gparts(struct part* const parts,
gparts[i + Ndm].mass = parts[i].mass;
/* Link the particles */
gparts[i + Ndm].part = &parts[i];
gparts[i + Ndm].id_or_neg_offset = -i;
parts[i].gpart = &gparts[i + Ndm];
}
}
......@@ -580,9 +578,8 @@ void collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
* gparts[i].part); */
/* And collect the DM ones */
if (gparts[i].id < 0LL) {
memcpy(&dmparts[count], &gparts[i], sizeof(struct gpart));
dmparts[count].id = -dmparts[count].id;
if (gparts[i].id_or_neg_offset > 0) {
dmparts[count] = gparts[i];
count++;
}
}
......
......@@ -59,8 +59,8 @@
*
* (Should be used for debugging only as it runs in O(N).)
*/
void printParticle(struct part *parts, struct xpart *xparts, long long int id,
size_t N) {
void printParticle(const struct part *parts, struct xpart *xparts,
long long int id, size_t N) {
int found = 0;
......@@ -82,24 +82,27 @@ void printParticle(struct part *parts, struct xpart *xparts, long long int id,
* the standard output.
*
* @param gparts The array of g-particles.
* @param parts The array of particles.
* @param id The id too look for.
* @param N The size of the array of g-particles.
*
* (Should be used for debugging only as it runs in O(N).)
*/
void printgParticle(struct gpart *gparts, long long int id, size_t N) {
void printgParticle(const struct gpart *gparts, const struct part *parts,
long long int id, size_t N) {
int found = 0;
/* Look for the particle. */
for (size_t i = 0; i < N; i++)
if (gparts[i].id == -id) {
printf("## gParticle[%zd] (DM) :\n id=%lld ", i, -gparts[i].id);
if (gparts[i].id_or_neg_offset == id) {
printf("## gParticle[%zd] (DM) :\n id=%lld", i, id);
gravity_debug_particle(&gparts[i]);
found = 1;
break;
} else if (gparts[i].id > 0 && gparts[i].part->id == id) {
printf("## gParticle[%zd] (hydro) :\n id=%lld ", i, gparts[i].id);
} else if (gparts[i].id_or_neg_offset < 0 &&
parts[-gparts[i].id_or_neg_offset].id == id) {
printf("## gParticle[%zd] (hydro) :\n id=%lld", i, id);
gravity_debug_particle(&gparts[i]);
found = 1;
break;
......@@ -114,9 +117,9 @@ void printgParticle(struct gpart *gparts, long long int id, size_t N) {
* @param p The particle to print
* @param xp The extended data ot the particle to print
*/
void printParticle_single(struct part *p, struct xpart *xp) {
void printParticle_single(const struct part *p, const struct xpart *xp) {
printf("## Particle: id=%lld ", p->id);
printf("## Particle: id=%lld", p->id);
hydro_debug_particle(p, xp);
printf("\n");
}
......@@ -128,7 +131,7 @@ void printParticle_single(struct part *p, struct xpart *xp) {
*/
void printgParticle_single(struct gpart *gp) {
printf("## g-Particle: id=%lld ", gp->id);
printf("## g-Particle: id=%lld ", gp->id_or_neg_offset);
gravity_debug_particle(gp);
printf("\n");
}
......
......@@ -19,14 +19,15 @@
#ifndef SWIFT_DEBUG_H
#define SWIFT_DEBUG_H
struct part;
struct gpart;
struct xpart;
/* Includes. */
#include "cell.h"
#include "part.h"
void printParticle(struct part *parts, struct xpart *xparts, long long int id,
size_t N);
void printgParticle(struct gpart *parts, long long int id, size_t N);
void printParticle_single(struct part *p, struct xpart *xp);
void printParticle(const struct part *parts, struct xpart *xparts,
long long int id, size_t N);
void printgParticle(const struct gpart *gparts, const struct part *parts,
long long int id, size_t N);
void printParticle_single(const struct part *p, const struct xpart *xp);
void printgParticle_single(struct gpart *gp);
#ifdef HAVE_METIS
......
......@@ -331,11 +331,11 @@ void engine_redistribute(struct engine *e) {
}
#ifdef SWIFT_DEBUG_CHECKS
if (s->parts[k].gpart->id < 0)
if (s->parts[k].gpart->id_or_neg_offset >= 0)
error("Trying to link a partnerless gpart !");
#endif
s->parts[k].gpart->id = count_this_dest;
s->parts[k].gpart->id_or_neg_offset = -count_this_dest;
count_this_dest++;
}
}
......@@ -519,13 +519,14 @@ void engine_redistribute(struct engine *e) {
for (size_t k = offset_gparts; k < offset_gparts + count_gparts; ++k) {
/* Does this gpart have a partner ? */
if (gparts_new[k].id >= 0) {
if (gparts_new[k].id_or_neg_offset <= 0) {
const size_t partner_index = offset_parts + gparts_new[k].id;
const ptrdiff_t partner_index =
offset_parts - gparts_new[k].id_or_neg_offset;
/* Re-link */
gparts_new[k].part = &parts_new[partner_index];
gparts_new[k].part->gpart = &gparts_new[k];
gparts_new[k].id_or_neg_offset = -partner_index;
parts_new[partner_index].gpart = &gparts_new[k];
}
}
......@@ -546,22 +547,22 @@ void engine_redistribute(struct engine *e) {
/* Verify that the links are correct */
for (size_t k = 0; k < nr_gparts; ++k) {
if (gparts_new[k].id > 0) {
if (gparts_new[k].id_or_neg_offset <= 0) {
if (gparts_new[k].part->gpart != &gparts_new[k])
error("Linking problem !");
struct part *part = &parts_new[-gparts_new[k].id_or_neg_offset];
if (gparts_new[k].x[0] != gparts_new[k].part->x[0] ||
gparts_new[k].x[1] != gparts_new[k].part->x[1] ||
gparts_new[k].x[2] != gparts_new[k].part->x[2])
if (part->gpart != &gparts_new[k]) error("Linking problem !");
if (gparts_new[k].x[0] != part->x[0] ||
gparts_new[k].x[1] != part->x[1] || gparts_new[k].x[2] != part->x[2])
error("Linked particles are not at the same position !");
}
}
for (size_t k = 0; k < nr_parts; ++k) {
if (parts_new[k].gpart != NULL) {
if (parts_new[k].gpart->part != &parts_new[k]) error("Linking problem !");
if (parts_new[k].gpart != NULL &&
parts_new[k].gpart->id_or_neg_offset != -k) {
error("Linking problem !");
}
}
#endif
......@@ -725,7 +726,6 @@ void engine_addtasks_send(struct engine *e, struct cell *ci, struct cell *cj) {
* @param t_xv The recv_xv #task, if it has already been created.
* @param t_rho The recv_rho #task, if it has already been created.
*/
void engine_addtasks_recv(struct engine *e, struct cell *c, struct task *t_xv,
struct task *t_rho) {
......@@ -781,7 +781,8 @@ void engine_exchange_cells(struct engine *e) {
MPI_Status status;
ticks tic = getticks();
/* Run through the cells and get the size of the ones that will be sent */
/* Run through the cells and get the size of the ones that will be sent off.
*/
int count_out = 0;
for (int k = 0; k < nr_cells; k++) {
offset[k] = count_out;
......@@ -946,7 +947,8 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
/* Re-link the associated gpart with the buffer offset of the part. */
if (s->parts[offset_parts + k].gpart != NULL) {
s->parts[offset_parts + k].gpart->id = e->proxies[pid].nr_parts_out;
s->parts[offset_parts + k].gpart->id_or_neg_offset =
-e->proxies[pid].nr_parts_out;
}
/* Load the part and xpart into the proxy. */
......@@ -962,7 +964,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
error(
"Do not have a proxy for the requested nodeID %i for part with "
"id=%lli, x=[%e,%e,%e].",
node_id, s->gparts[offset_parts + k].id,
node_id, s->gparts[offset_parts + k].id_or_neg_offset,
s->gparts[offset_gparts + k].x[0], s->gparts[offset_parts + k].x[1],
s->gparts[offset_gparts + k].x[2]);
proxy_gparts_load(&e->proxies[pid], &s->gparts[offset_gparts + k], 1);
......@@ -1022,7 +1024,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
s->xparts = xparts_new;
for (size_t k = 0; k < offset_parts; k++) {
if (s->parts[k].gpart != NULL) {
s->parts[k].gpart->part = &s->parts[k];
s->parts[k].gpart->id_or_neg_offset = -k;
}
}
}
......@@ -1037,8 +1039,8 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
free(s->gparts);
s->gparts = gparts_new;
for (size_t k = 0; k < offset_gparts; k++) {
if (s->gparts[k].id > 0) {
s->gparts[k].part->gpart = &s->gparts[k];
if (s->gparts[k].id_or_neg_offset < 0) {
s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
}
}
}
......@@ -1111,9 +1113,10 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
/* Re-link the gparts. */
for (int k = 0; k < p->nr_gparts_in; k++) {
struct gpart *gp = &s->gparts[offset_gparts + count_gparts + k];
if (gp->id >= 0) {
struct part *p = &s->parts[offset_parts + count_parts + gp->id];
gp->part = p;
if (gp->id_or_neg_offset <= 0) {
struct part *p =
&s->parts[offset_gparts + count_parts - gp->id_or_neg_offset];
gp->id_or_neg_offset = s->parts - p;
p->gpart = gp;
}
}
......@@ -1762,7 +1765,6 @@ void engine_print_task_counts(struct engine *e) {
*
* @param e The #engine.
*/
void engine_rebuild(struct engine *e) {
const ticks tic = getticks();
......@@ -2494,8 +2496,7 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
s->xparts = xparts_new;
/* Re-link the gparts. */
for (size_t k = 0; k < s->nr_parts; k++)
if (s->parts[k].gpart != NULL) s->parts[k].gpart->part = &s->parts[k];
part_relink_gparts(s->parts, s->nr_parts, 0);
/* Re-allocate the local gparts. */
if (e->verbose)
......@@ -2511,31 +2512,28 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
s->gparts = gparts_new;
/* Re-link the parts. */
for (size_t k = 0; k < s->nr_gparts; k++)
if (s->gparts[k].id > 0) s->gparts[k].part->gpart = &s->gparts[k];
part_relink_parts(s->gparts, s->nr_gparts, s->parts);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the links are correct */
for (size_t k = 0; k < s->nr_gparts; ++k) {
if (s->gparts[k].id > 0) {
if (s->gparts[k].id_or_neg_offset <= 0) {
struct part *part = &s->parts[-s->gparts[k].id_or_neg_offset];
if (s->gparts[k].part->gpart != &s->gparts[k]) error("Linking problem !");
if (part->gpart != &s->gparts[k]) error("Linking problem !");
if (s->gparts[k].x[0] != s->gparts[k].part->x[0] ||
s->gparts[k].x[1] != s->gparts[k].part->x[1] ||
s->gparts[k].x[2] != s->gparts[k].part->x[2])
if (s->gparts[k].x[0] != part->x[0] || s->gparts[k].x[1] != part->x[1] ||
s->gparts[k].x[2] != part->x[2])
error("Linked particles are not at the same position !");
}
}
for (size_t k = 0; k < s->nr_parts; ++k) {
if (s->parts[k].gpart != NULL) {
if (s->parts[k].gpart->part != &s->parts[k]) error("Linking problem !");
}
if (s->parts[k].gpart != NULL && s->parts[k].gpart->id_or_neg_offset != -k)
error("Linking problem !");
}
#endif
#else
......
......@@ -18,7 +18,7 @@
******************************************************************************/
__attribute__((always_inline)) INLINE static void gravity_debug_particle(
struct gpart* p) {
const struct gpart* p) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
......
......@@ -39,8 +39,8 @@ __attribute__((always_inline)) INLINE static void darkmatter_read_particles(
COMPULSORY);
readArray(h_grp, "Velocities", FLOAT, N, 3, gparts, N_total, offset, v_full,
COMPULSORY);
readArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, gparts, N_total, offset, id,
COMPULSORY);
readArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, gparts, N_total, offset,
id_or_neg_offset, COMPULSORY);
}
/**
......@@ -75,6 +75,6 @@ __attribute__((always_inline)) INLINE static void darkmatter_write_particles(
Ndm, 3, gparts, Ndm_total, mpi_rank, offset, v_full, us,
UNIT_CONV_SPEED);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "ParticleIDs",
ULONGLONG, Ndm, 1, gparts, Ndm_total, mpi_rank, offset, id, us,
UNIT_CONV_NO_UNITS);
ULONGLONG, Ndm, 1, gparts, Ndm_total, mpi_rank, offset,
id_or_neg_offset, us, UNIT_CONV_NO_UNITS);
}
......@@ -47,14 +47,8 @@ struct gpart {
/* float tx; */
/* float tv; */
/* Anonymous union for id/part. */
union {
/* Particle ID. */
long long id;
/* Pointer to corresponding SPH part. */
struct part* part;
};
/* Particle ID. If negative, it is the negative offset of the #part with
which this gpart is linked. */
long long id_or_neg_offset;
} __attribute__((aligned(gpart_align)));
......@@ -18,7 +18,7 @@
******************************************************************************/
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
struct part* p, struct xpart* xp) {
const struct part* p, const struct xpart* xp) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
......
......@@ -18,7 +18,7 @@
******************************************************************************/
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
struct part* p, struct xpart* xp) {
const struct part* p, const struct xpart* xp) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
......
......@@ -18,7 +18,7 @@
******************************************************************************/
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
struct part* p, struct xpart* xp) {
const struct part* p, const struct xpart* xp) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
......
......@@ -29,6 +29,36 @@
#include "error.h"
#include "part.h"
/**
* @brief Re-link the #gpart%s associated with the list of #part%s.
*
* @param parts The list of #part.
* @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) {
for (size_t k = 0; k < N; k++) {
if (parts[k].gpart) {
parts[k].gpart->id_or_neg_offset = -(k + offset);
}
}
}
/**
* @brief Re-link the #gpart%s associated with the list of #part%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.
*/
void part_relink_parts(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) {
parts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
}
}
}
#ifdef WITH_MPI
/* MPI data type for the particle transfers */
MPI_Datatype part_mpi_type;
......
......@@ -22,6 +22,9 @@
/* Config parameters. */
#include "../config.h"
/* Standard headers. */
#include <stddef.h>
/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
......@@ -49,6 +52,8 @@
/* Import the right gravity particle definition */
#include "./gravity/Default/gravity_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);
#ifdef WITH_MPI
/* MPI data type for the particle transfers */
extern MPI_Datatype part_mpi_type;
......
......@@ -760,7 +760,7 @@ void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
struct gpart *const gp = &gparts[k];
/* If the g-particle has no counterpart */
if (gp->id < 0) {
if (gp->id_or_neg_offset > 0) {
/* First, finish the force calculation */
gravity_end_force(gp);
......@@ -872,7 +872,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
struct gpart *const gp = &gparts[k];
/* If the g-particle has no counterpart and needs to be kicked */
if (gp->id < 0) {
if (gp->id_or_neg_offset > 0) {
if (gp->ti_end <= ti_current) {
......
......@@ -442,10 +442,10 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
s->parts[k] = s->parts[nr_parts];
s->parts[nr_parts] = tp;
if (s->parts[k].gpart != NULL) {
s->parts[k].gpart->part = &s->parts[k];
s->parts[k].gpart->id_or_neg_offset = -k;
}
if (s->parts[nr_parts].gpart != NULL) {
s->parts[nr_parts].gpart->part = &s->parts[nr_parts];
s->parts[nr_parts].gpart->id_or_neg_offset = -nr_parts;
}
const struct xpart txp = s->xparts[k];
s->xparts[k] = s->xparts[nr_parts];
......@@ -481,11 +481,12 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
const struct gpart tp = s->gparts[k];
s->gparts[k] = s->gparts[nr_gparts];
s->gparts[nr_gparts] = tp;
if (s->gparts[k].id > 0) {
s->gparts[k].part->gpart = &s->gparts[k];
if (s->gparts[k].id_or_neg_offset <= 0) {
s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
}
if (s->gparts[nr_gparts].id > 0) {
s->gparts[nr_gparts].part->gpart = &s->gparts[nr_gparts];
if (s->gparts[nr_gparts].id_or_neg_offset <= 0) {
s->parts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
&s->gparts[nr_gparts];
}
const int t = gind[k];
gind[k] = gind[nr_gparts];
......@@ -551,8 +552,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
space_parts_sort(s, ind, nr_parts, 0, s->nr_cells - 1, verbose);
/* Re-link the gparts. */
for (size_t k = 0; k < nr_parts; k++)
if (s->parts[k].gpart != NULL) s->parts[k].gpart->part = &s->parts[k];
part_relink_gparts(s->parts, nr_parts, 0);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify space_sort_struct. */
......@@ -600,8 +600,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
space_gparts_sort(s, gind, nr_gparts, 0, s->nr_cells - 1, verbose);
/* Re-link the parts. */
for (int k = 0; k < nr_gparts; k++)
if (s->gparts[k].id > 0) s->gparts[k].part->gpart = &s->gparts[k];
part_relink_parts(s->gparts, nr_gparts, s->parts);
/* We no longer need the indices as of here. */
free(gind);
......@@ -610,21 +609,22 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
/* Verify that the links are correct */
for (size_t k = 0; k < nr_gparts; ++k) {
if (s->gparts[k].id > 0) {
if (s->gparts[k].id_or_neg_offset < 0) {
if (s->gparts[k].part->gpart != &s->gparts[k]) error("Linking problem !");
const struct part *part = &s->parts[-s->gparts[k].id_or_neg_offset];
if (s->gparts[k].x[0] != s->gparts[k].part->x[0] ||
s->gparts[k].x[1] != s->gparts[k].part->x[1] ||
s->gparts[k].x[2] != s->gparts[k].part->x[2])
if (part->gpart != &s->gparts[k]) error("Linking problem !");
if (s->gparts[k].x[0] != part->x[0] || s->gparts[k].x[1] != part->x[1] ||
s->gparts[k].x[2] != part->x[2])
error("Linked particles are not at the same position !");
}
}
for (size_t k = 0; k < nr_parts; ++k) {
if (s->parts[k].gpart != NULL) {