diff --git a/examples/main.c b/examples/main.c
index 389fc52ee19ffbfe818e678f9ceee78752e124b7..97ad5dd4638be0b69e95661acfef6af6ed11c655 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -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;
}
diff --git a/src/cell.c b/src/cell.c
index 1531c33e827c8514b14ef7d1d54e123bef1842ad..d1952ff7ea86a541ea766520ff9b9b41bacb8d49 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -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);
}
/**
diff --git a/src/cell.h b/src/cell.h
index 31efb179957f462bc12c2681973a183f8e7aa269..89107c3d1b4a8e3b98112cdb48de3986bcc4ce7e 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -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);
diff --git a/src/common_io.c b/src/common_io.c
index 03d5618cd0cafe5ff6067cf72107baacffb5d85a..971fe6b01c682c489d3444ec1d47d6a902250bb8 100644
--- a/src/common_io.c
+++ b/src/common_io.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++;
}
}
diff --git a/src/debug.c b/src/debug.c
index 77238ab135e3f27b8c2c43b0ec18e7745493b138..e7b60c52f92fdb64fdefc8ac5c5e353f1ab57e6e 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -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");
}
diff --git a/src/debug.h b/src/debug.h
index 13be15adb867e8bafe95e2900f68caaa36192510..367241201977d9b79a8c2913dbae5d08f1148529 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -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
diff --git a/src/engine.c b/src/engine.c
index 79bcd3ec2f84f91fb65108851911b3427b21045d..7f4361d07a4e98350c9687bfda7aef85208e7adb 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -332,11 +332,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++;
}
}
@@ -520,13 +520,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];
}
}
@@ -547,22 +548,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
@@ -947,7 +948,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. */
@@ -963,7 +965,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);
@@ -1023,7 +1025,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;
}
}
}
@@ -1038,8 +1040,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];
}
}
}
@@ -1112,9 +1114,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;
}
}
@@ -2686,8 +2689,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)
@@ -2703,31 +2705,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) {
- if (s->gparts[k].part->gpart != &s->gparts[k]) error("Linking problem !");
+ 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 < 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
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
index 21d775d703c8bd7000ad03386fe818124642a8f1..7cf375a1fdf7bccc4131dc415ab2d4acbbf2d3bc 100644
--- a/src/gravity/Default/gravity_debug.h
+++ b/src/gravity/Default/gravity_debug.h
@@ -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 "
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
index 129c4b39828ca73d2d80d79edbdaa8ec4d5a9e01..74f364dd97361f0513755bedec83fe7cb277c36b 100644
--- a/src/gravity/Default/gravity_io.h
+++ b/src/gravity/Default/gravity_io.h
@@ -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);
}
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index afa207f9e78a839a7b2bca2f328fcdde432c3cc7..3001dd0406d1c6fa1047437e00ac1a80eb535b12 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -48,14 +48,8 @@ struct gpart {
double mass_interacted;
- /* 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)));
diff --git a/src/hydro/Default/hydro_debug.h b/src/hydro/Default/hydro_debug.h
index 79ee392d46ca75a2c097bf045b2d82c9f3dc96c0..9b8136e8349398e4e9e804459a5de23d21043564 100644
--- a/src/hydro/Default/hydro_debug.h
+++ b/src/hydro/Default/hydro_debug.h
@@ -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 "
diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h
index 31e89d438bc96c0f0f2ba56249664d28036cb607..b67e79182ccaaee7c0421c57a91ec9fa2adae65c 100644
--- a/src/hydro/Gadget2/hydro_debug.h
+++ b/src/hydro/Gadget2/hydro_debug.h
@@ -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 "
diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h
index 127ba75e99418b6a5dc197a44ccdc77de3cdef15..1e0e3f1f111149a3babeef431893dff2cbe6bb3c 100644
--- a/src/hydro/Minimal/hydro_debug.h
+++ b/src/hydro/Minimal/hydro_debug.h
@@ -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], "
diff --git a/src/part.c b/src/part.c
index b89abdde40fe8c7a57d1e9ac9e18fece83ba1f21..b00eaccaae0e86f7c4e8019a307f0bf455687b7c 100644
--- a/src/part.c
+++ b/src/part.c
@@ -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;
diff --git a/src/part.h b/src/part.h
index e99be6e51a9bd74dd9eec8f590e80989e83ec2e1..efca7b6b5bef49f20df1e2c45b30f65ecbbf4960 100644
--- a/src/part.h
+++ b/src/part.h
@@ -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;
diff --git a/src/runner.c b/src/runner.c
index a263ebb03866d2e10c37b950f2733b9bf8581458..88df1824f832b2a1294e3e6b0e08f615120c58cd 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -769,7 +769,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, const_G);
@@ -882,7 +882,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) {
diff --git a/src/space.c b/src/space.c
index cd5cf1179cc9da2d0a89b18ce54927ef360e7870..12147cf298cf28ec5737340b72aca65b290b74c9 100644
--- a/src/space.c
+++ b/src/space.c
@@ -448,10 +448,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];
@@ -487,11 +487,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];
@@ -557,8 +558,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. */
@@ -606,8 +606,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);
@@ -616,21 +615,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) {
-
- 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
@@ -1274,7 +1274,7 @@ void space_do_split(struct space *s, struct cell *c) {
}
/* Split the cell data. */
- cell_split(c);
+ cell_split(c, c->parts - s->parts);
/* Remove any progeny with zero parts. */
for (int k = 0; k < 8; k++)
diff --git a/src/tools.c b/src/tools.c
index a60c785104a0bb1bb5e508954fb4ee68f5106668..568adcdb680f88dcdd0c456c0ef6f63a6b70f38b 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -181,16 +181,19 @@ void pairs_single_density(double *dim, long long int pid,
}
void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
+
float r2, hi, hj, hig2, hjg2, dx[3];
struct part *pi, *pj;
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->count; ++i) {
+
pi = &ci->parts[i];
hi = pi->h;
hig2 = hi * hi * kernel_gamma2;
for (int j = 0; j < cj->count; ++j) {
+
pj = &cj->parts[j];
/* Pairwise distance */
@@ -202,6 +205,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
/* Hit or miss? */
if (r2 < hig2) {
+
/* Interact */
runner_iact_nonsym_density(r2, dx, hi, pj->h, pi, pj);
}
@@ -210,11 +214,13 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
/* Reverse double-for loop and checks every interaction */
for (int j = 0; j < cj->count; ++j) {
+
pj = &cj->parts[j];
hj = pj->h;
hjg2 = hj * hj * kernel_gamma2;
for (int i = 0; i < ci->count; ++i) {
+
pi = &ci->parts[i];
/* Pairwise distance */
@@ -226,6 +232,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
/* Hit or miss? */
if (r2 < hjg2) {
+
/* Interact */
runner_iact_nonsym_density(r2, dx, hj, pi->h, pj, pi);
}
@@ -239,11 +246,13 @@ void self_all_density(struct runner *r, struct cell *ci) {
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->count; ++i) {
+
pi = &ci->parts[i];
hi = pi->h;
hig2 = hi * hi * kernel_gamma2;
for (int j = i + 1; j < ci->count; ++j) {
+
pj = &ci->parts[j];
hj = pj->h;
hjg2 = hj * hj * kernel_gamma2;
@@ -259,12 +268,14 @@ void self_all_density(struct runner *r, struct cell *ci) {
/* Hit or miss? */
if (r2 < hig2) {
+
/* Interact */
runner_iact_nonsym_density(r2, dxi, hi, hj, pi, pj);
}
/* Hit or miss? */
if (r2 < hjg2) {
+
dxi[0] = -dxi[0];
dxi[1] = -dxi[1];
dxi[2] = -dxi[2];
@@ -277,7 +288,9 @@ void self_all_density(struct runner *r, struct cell *ci) {
}
void pairs_single_grav(double *dim, long long int pid,
- struct gpart *restrict parts, int N, int periodic) {
+ struct gpart *restrict gparts, const struct part *parts,
+ int N, int periodic) {
+
int i, k;
// int mj, mk;
// double maxratio = 1.0;
@@ -288,18 +301,20 @@ void pairs_single_grav(double *dim, long long int pid,
/* Find "our" part. */
for (k = 0; k < N; k++)
- if ((parts[k].id > 0 && parts[k].part->id == pid) || parts[k].id == -pid)
+ if ((gparts[k].id_or_neg_offset < 0 &&
+ parts[-gparts[k].id_or_neg_offset].id == pid) ||
+ gparts[k].id_or_neg_offset == pid)
break;
if (k == N) error("Part not found.");
- pi = parts[k];
+ pi = gparts[k];
pi.a_grav[0] = 0.0f;
pi.a_grav[1] = 0.0f;
pi.a_grav[2] = 0.0f;
/* Loop over all particle pairs. */
for (k = 0; k < N; k++) {
- if (parts[k].id == pi.id) continue;
- pj = parts[k];
+ if (gparts[k].id_or_neg_offset == pi.id_or_neg_offset) continue;
+ pj = gparts[k];
for (i = 0; i < 3; i++) {
dx[i] = pi.x[i] - pj.x[i];
if (periodic) {
@@ -326,7 +341,8 @@ void pairs_single_grav(double *dim, long long int pid,
/* Dump the result. */
message(
"acceleration on gpart %lli is a=[ %e %e %e ], |a|=[ %.2e %.2e %.2e ].\n",
- pi.part->id, a[0], a[1], a[2], aabs[0], aabs[1], aabs[2]);
+ parts[-pi.id_or_neg_offset].id, a[0], a[1], a[2], aabs[0], aabs[1],
+ aabs[2]);
}
/**
@@ -335,6 +351,7 @@ void pairs_single_grav(double *dim, long long int pid,
* @param N number of intervals in [0,1].
*/
void density_dump(int N) {
+
int k;
float r2[4] = {0.0f, 0.0f, 0.0f, 0.0f}, hi[4], hj[4];
struct part /**pi[4], *pj[4],*/ Pi[4], Pj[4];
diff --git a/src/tools.h b/src/tools.h
index a2320b7690499c06edee96bdb008cee7b707c19e..e2a9497ef747fa96cee2813b5950b2300f1b0dd9 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -30,7 +30,8 @@
void factor(int value, int *f1, int *f2);
void density_dump(int N);
void pairs_single_grav(double *dim, long long int pid,
- struct gpart *restrict parts, int N, int periodic);
+ struct gpart *restrict gparts, const struct part *parts,
+ int N, int periodic);
void pairs_single_density(double *dim, long long int pid,
struct part *restrict parts, int N, int periodic);