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

src/cell.c

@@ -451,7 +451,9 @@ 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];
+    if (parts[k].gpart != NULL) {
+      parts[k].gpart->id_or_neg_offset = -k;
+    }
 
   /* Verify that _all_ the parts have been assigned to a cell. */
   /* for ( k = 1 ; k < 8 ; k++ )
@@ -545,7 +547,9 @@ 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];
+    if (gparts[k].id_or_neg_offset < 0) {
+      parts[gparts[k].id_or_neg_offset].gpart = &gparts[k];
+    }
 }
 
 /**

src/common_io.c

@@ -489,9 +489,8 @@ void prepare_dm_gparts(struct gpart* gparts, size_t Ndm) {
   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");
-
-    gparts[i].id = -gparts[i].id;
+    if (gparts[i].id_or_neg_offset <= 0)
+      error("0 or negative ID for DM particle");
   }
 }
 
@@ -524,7 +523,7 @@ void duplicate_hydro_gparts(struct part* parts, struct gpart* gparts,
     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];
   }
 }
@@ -546,9 +545,8 @@ void collect_dm_gparts(struct gpart* gparts, size_t Ntot, struct gpart* dmparts,
   for (size_t i = 0; i < Ntot; ++i) {
 
     /* And collect the DM ones */
-    if (gparts[i].id < 0) {
+    if (gparts[i].id_or_neg_offset > 0) {
       memcpy(&dmparts[count], &gparts[i], sizeof(struct gpart));
-      dmparts[count].id = -dmparts[count].id;
       count++;
     }
   }

src/debug.c

@@ -52,8 +52,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;
 
@@ -69,19 +69,21 @@ void printParticle(struct part *parts, struct xpart *xparts, long long int id,
   if (!found) printf("## Particles[???] id=%lld not found\n", id);
 }
 
-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;
@@ -98,7 +100,7 @@ void printgParticle(struct gpart *gparts, long long int id, size_t N) {
  *
  */
 
-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);
   hydro_debug_particle(p, xp);

src/debug.h

@@ -23,10 +23,11 @@
 #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);
 
 #ifdef HAVE_METIS
 #include "metis.h"

src/gravity/Default/gravity_debug.h

@@ -18,7 +18,7 @@
  ******************************************************************************/
 
 __attribute__((always_inline))
-    INLINE static void gravity_debug_particle(struct gpart* p) {
+    INLINE static void gravity_debug_particle(const struct gpart* p) {
   printf(
       "x=[%.3e,%.3e,%.3e], "
       "v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "

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);
 }
 
 /**
@@ -71,5 +71,6 @@ __attribute__((always_inline)) INLINE static void darkmatter_write_particles(
   writeArray(h_grp, fileName, xmfFile, "Velocities", FLOAT, Ndm, 3, gparts,
              Ndm_total, mpi_rank, offset, v_full, us, UNIT_CONV_SPEED);
   writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, Ndm, 1, gparts,
-             Ndm_total, mpi_rank, offset, id, us, UNIT_CONV_NO_UNITS);
+             Ndm_total, mpi_rank, offset, id_or_neg_offset, us,
+             UNIT_CONV_NO_UNITS);
 }

src/gravity/Default/gravity_part.h

@@ -40,14 +40,8 @@ struct gpart {
   /* Particle time of end of time-step. */
   int ti_end;
 
-  /* 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(part_align)));

src/hydro/Gadget2/hydro_debug.h

@@ -17,8 +17,8 @@
  *
  ******************************************************************************/
 
-__attribute__((always_inline))
-    INLINE static void hydro_debug_particle(struct part* p, struct xpart* xp) {
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+    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 "

src/space.c

@@ -403,8 +403,10 @@ 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 (int k = 0; k < nr_parts; k++)
-    if (s->parts[k].gpart != NULL) s->parts[k].gpart->part = &s->parts[k];
+  for (size_t k = 0; k < nr_parts; k++)
+    if (s->parts[k].gpart != NULL) {
+      s->parts[k].gpart->id_or_neg_offset = -k;
+    }
 
   /* Verify space_sort_struct. */
   /* for ( k = 1 ; k < nr_parts ; k++ ) {
@@ -491,7 +493,9 @@ void space_rebuild(struct space *s, double cell_max, int 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];
+    if (s->gparts[k].id_or_neg_offset < 0) {
+      s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
+    }
 
   /* We no longer need the indices as of here. */
   free(gind);

src/tools.c

@@ -237,7 +237,8 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
 }
 
 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;
@@ -249,18 +250,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[0] = 0.0f;
   pi.a[1] = 0.0f;
   pi.a[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) {
@@ -287,7 +290,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]);
 }
 
 /**

src/tools.h

@@ -28,7 +28,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);