Add test27cellsStars

.gitignore

@@ -63,12 +63,18 @@ tests/brute_force_perturbed.dat
 tests/swift_dopair_perturbed.dat
 tests/test27cells
 tests/test27cells_subset
+tests/test27cellsStars
+tests/test27cellsStars_subset
 tests/testPeriodicBC
 tests/test125cells
 tests/brute_force_27_standard.dat
 tests/swift_dopair_27_standard.dat
 tests/brute_force_27_perturbed.dat
 tests/swift_dopair_27_perturbed.dat
+tests/star_brute_force_27_standard.dat
+tests/swift_star_dopair_27_standard.dat
+tests/star_brute_force_27_perturbed.dat
+tests/swift_star_dopair_27_perturbed.dat
 tests/brute_force_125_standard.dat
 tests/swift_dopair_125_standard.dat
 tests/brute_force_125_perturbed.dat
@@ -106,6 +112,8 @@ tests/testPeriodicBC.sh
 tests/testPeriodicBCPerturbed.sh
 tests/test27cells.sh
 tests/test27cellsPerturbed.sh
+tests/test27cellsStars.sh
+tests/test27cellsStarsPerturbed.sh
 tests/test125cells.sh
 tests/test125cellsPerturbed.sh
 tests/testParser.sh

configure.ac

@@ -1427,6 +1427,8 @@ AC_CONFIG_FILES([tests/testReading.sh], [chmod +x tests/testReading.sh])
 AC_CONFIG_FILES([tests/testActivePair.sh], [chmod +x tests/testActivePair.sh])
 AC_CONFIG_FILES([tests/test27cells.sh], [chmod +x tests/test27cells.sh])
 AC_CONFIG_FILES([tests/test27cellsPerturbed.sh], [chmod +x tests/test27cellsPerturbed.sh])
+AC_CONFIG_FILES([tests/test27cellsStars.sh], [chmod +x tests/test27cellsStars.sh])
+AC_CONFIG_FILES([tests/test27cellsStarsPerturbed.sh], [chmod +x tests/test27cellsStarsPerturbed.sh])
 AC_CONFIG_FILES([tests/test125cells.sh], [chmod +x tests/test125cells.sh])
 AC_CONFIG_FILES([tests/test125cellsPerturbed.sh], [chmod +x tests/test125cellsPerturbed.sh])
 AC_CONFIG_FILES([tests/testPeriodicBC.sh], [chmod +x tests/testPeriodicBC.sh])

src/cell.c

@@ -61,6 +61,7 @@
 #include "scheduler.h"
 #include "space.h"
 #include "space_getsid.h"
+#include "stars.h"
 #include "timers.h"
 #include "tools.h"
 
@@ -2874,6 +2875,8 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
 
       /* Drift... */
       drift_spart(sp, dt_drift, ti_old_gpart, ti_current);
+      if (spart_is_active(sp, e))
+	star_init_spart(sp);
 
       /* Note: no need to compute dx_max as all spart have a gpart */
     }

src/runner.c

@@ -344,12 +344,6 @@ void runner_do_star_ghost(struct runner *r, struct cell *c, int timer) {
     free(sid);
   }
 
-  /* Reset the data for next loop */
-  for (int k = 0; k < c->scount; k++)
-    if (spart_is_active(&sparts[k], e)) {
-      star_init_spart(&sparts[k]);
-    }
-
   if (timer) TIMER_TOC(timer_do_star_ghost);
 }
 

src/runner_doiact_star.h

@@ -78,7 +78,7 @@ void runner_doself_star_density(struct runner *r, struct cell *c, int timer) {
         error("Particle pj not drifted to current time");
 #endif
 
-      if (r2 < hig2) {
+      if (r2 > 0.f && r2 < hig2) {
 	runner_iact_nonsym_star_density(r2, dx, hi, hj, si, pj, a, H);
       }
     } /* loop over the parts in ci. */
@@ -89,7 +89,7 @@ void runner_doself_star_density(struct runner *r, struct cell *c, int timer) {
 }
 
 /**
- * @brief Calculate the number density of #part around the #spart
+ * @brief Calculate the number density of cj #part around the ci #spart
  *
  * @param r runner task
  * @param c cell
@@ -968,9 +968,6 @@ void runner_doself_branch_star_density(struct runner *r, struct cell *c) {
   if (c->h_max_old * kernel_gamma > c->dmin)
     error("Cell smaller than smoothing length");
 
-  /* /\* Check that cells are drifted. *\/ */
-  /* if (!cell_are_part_drifted(c, e)) error("Interacting undrifted cell."); */
-
   runner_doself_star_density(r, c, 1);
 }
 
@@ -1281,9 +1278,9 @@ void runner_dosub_pair_star_density(struct runner *r, struct cell *ci, struct ce
   /* Otherwise, compute the pair directly. */
   else if (cell_is_active_star(ci, e) || cell_is_active_star(cj, e)) {
 
-    /* /\* Make sure both cells are drifted to the current timestep. *\/ */
-    /* if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e)) */
-    /*   error("Interacting undrifted cells."); */
+    /* Make sure both cells are drifted to the current timestep. */
+    if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
+      error("Interacting undrifted cells.");
 
     /* Do any of the cells need to be sorted first? */
     if (!(ci->sorted & (1 << sid)) ||
@@ -1329,10 +1326,7 @@ void runner_dosub_self_star_density(struct runner *r, struct cell *ci, int getti
   }
 
   /* Otherwise, compute self-interaction. */
-  else {
-
-    /* /\* Check that cells are drifted. *\/ */
-    /* if (!cell_are_part_drifted(c, e)) error("Interacting undrifted cell."); */
+ else {
 
     runner_doself_branch_star_density(r, ci);
   }

src/stars/Default/star.h

@@ -106,10 +106,6 @@ __attribute__((always_inline)) INLINE static void star_end_density(
   const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
   const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
 
-  /* Final operation on the density (add self-contribution). */
-  sp->wcount += kernel_root;
-  sp->wcount_dh -= hydro_dimension * kernel_root;
-
   /* Finish the calculation by inserting the missing h-factors */
   sp->wcount *= h_inv_dim;
   sp->wcount_dh *= h_inv_dim_plus_one;

src/tools.c

@@ -45,6 +45,7 @@
 #include "part.h"
 #include "periodic.h"
 #include "runner.h"
+#include "stars.h"
 
 /**
  *  Factorize a given integer, attempts to keep larger pair of factors.
@@ -335,6 +336,78 @@ void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj) {
   }
 }
 
+void pairs_all_star_density(struct runner *r, struct cell *ci, struct cell *cj) {
+
+  float r2, dx[3];
+  const double dim[3] = {r->e->s->dim[0], r->e->s->dim[1], r->e->s->dim[2]};
+  const struct engine *e = r->e;
+  const struct cosmology *cosmo = e->cosmology;
+  const float a = cosmo->a;
+  const float H = cosmo->H;
+
+  /* Implements a double-for loop and checks every interaction */
+  for (int i = 0; i < ci->scount; ++i) {
+    struct spart *spi = &ci->sparts[i];
+
+    float hi = spi->h;
+    float hig2 = hi * hi * kernel_gamma2;
+
+    /* Skip inactive particles. */
+    if (!spart_is_active(spi, e)) continue;
+
+    for (int j = 0; j < cj->count; ++j) {
+
+      struct part *pj = &cj->parts[j];
+
+      /* Pairwise distance */
+      r2 = 0.0f;
+      for (int k = 0; k < 3; k++) {
+        dx[k] = spi->x[k] - pj->x[k];
+        dx[k] = nearest(dx[k], dim[k]);
+        r2 += dx[k] * dx[k];
+      }
+
+      /* Hit or miss? */
+      if (r2 < hig2) {
+
+        /* Interact */
+        runner_iact_nonsym_star_density(r2, dx, hi, pj->h, spi, pj, a, H);
+      }
+    }
+  }
+
+  /* Reverse double-for loop and checks every interaction */
+  for (int j = 0; j < cj->scount; ++j) {
+
+    struct spart *spj = &cj->sparts[j];
+    float hj = spj->h;
+    float hjg2 = hj * hj * kernel_gamma2;
+
+    /* Skip inactive particles. */
+    if (!spart_is_active(spj, e)) continue;
+
+    for (int i = 0; i < ci->count; ++i) {
+
+      struct part *pi = &ci->parts[i];
+
+      /* Pairwise distance */
+      r2 = 0.0f;
+      for (int k = 0; k < 3; k++) {
+        dx[k] = spj->x[k] - pi->x[k];
+        dx[k] = nearest(dx[k], dim[k]);
+        r2 += dx[k] * dx[k];
+      }
+
+      /* Hit or miss? */
+      if (r2 < hjg2) {
+
+        /* Interact */
+        runner_iact_nonsym_star_density(r2, dx, hj, pi->h, spj, pi, a, H);
+      }
+    }
+  }
+}
+
 void self_all_density(struct runner *r, struct cell *ci) {
   float r2, hi, hj, hig2, hjg2, dxi[3];  //, dxj[3];
   struct part *pi, *pj;
@@ -428,6 +501,48 @@ void self_all_force(struct runner *r, struct cell *ci) {
   }
 }
 
+void self_all_star_density(struct runner *r, struct cell *ci) {
+  float r2, hi, hj, hig2, dxi[3];  //, dxj[3];
+  struct spart *spi;
+  struct part *pj;
+  const struct engine *e = r->e;
+  const struct cosmology *cosmo = e->cosmology;
+  const float a = cosmo->a;
+  const float H = cosmo->H;
+
+  /* Implements a double-for loop and checks every interaction */
+  for (int i = 0; i < ci->scount; ++i) {
+
+    spi = &ci->sparts[i];
+    hi = spi->h;
+    hig2 = hi * hi * kernel_gamma2;
+
+    if (!spart_is_active(spi, e))
+      continue;
+
+    for (int j = 0; j < ci->count; ++j) {
+
+      pj = &ci->parts[j];
+      hj = pj->h;
+
+      /* Pairwise distance */
+      r2 = 0.0f;
+      for (int k = 0; k < 3; k++) {
+        dxi[k] = spi->x[k] - pj->x[k];
+        r2 += dxi[k] * dxi[k];
+      }
+
+      /* Hit or miss? */
+      if (r2 > 0.f && r2 < hig2) {
+
+        /* Interact */
+        runner_iact_nonsym_star_density(r2, dxi, hi, hj, spi, pj, a, H);
+      }
+
+    }
+  }
+}
+
 /**
  * @brief Compute the force on a single particle brute-force.
  */
@@ -544,6 +659,23 @@ void shuffle_particles(struct part *parts, const int count) {
   }
 }
 
+/**
+ * @brief Randomly shuffle an array of sparticles.
+ */
+void shuffle_sparticles(struct spart *sparts, const int scount) {
+  if (scount > 1) {
+    for (int i = 0; i < scount - 1; i++) {
+      int j = i + random_uniform(0., (double)(scount - 1 - i));
+
+      struct spart sparticle = sparts[j];
+
+      sparts[j] = sparts[i];
+
+      sparts[i] = sparticle;
+    }
+  }
+}
+
 /**
  * @brief Compares two values based on their relative difference: |a - b|/|a +
  * b|

src/tools.h

@@ -40,11 +40,14 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj);
 void self_all_density(struct runner *r, struct cell *ci);
 void pairs_all_force(struct runner *r, struct cell *ci, struct cell *cj);
 void self_all_force(struct runner *r, struct cell *ci);
+void pairs_all_star_density(struct runner *r, struct cell *ci, struct cell *cj);
+void self_all_star_density(struct runner *r, struct cell *ci);
 
 void pairs_n2(double *dim, struct part *restrict parts, int N, int periodic);
 
 double random_uniform(double a, double b);
 void shuffle_particles(struct part *parts, const int count);
+void shuffle_sparticles(struct spart *sparts, const int scount);
 void gravity_n2(struct gpart *gparts, const int gcount,
                 const struct phys_const *constants,
                 const struct gravity_props *gravity_properties, float rlr);

tests/Makefile.am

@@ -28,7 +28,8 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
         testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
 	testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
 	testPotentialPair testEOS testUtilities testSelectOutput.sh \
-	testCbrt testCosmology testOutputList testFormat.sh
+	testCbrt testCosmology testOutputList testFormat.sh \
+	testCbrt testCosmology test27cellsStars.sh test27cellsStarsPerturbed.sh
 
 # List of test programs to compile
 check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
@@ -39,7 +40,8 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
                  testRiemannHLLC testMatrixInversion testDump testLogger \
 		 testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
 		 testGravityDerivatives testPotentialSelf testPotentialPair testEOS testUtilities \
-		 testSelectOutput testCbrt testCosmology testOutputList
+		 testSelectOutput testCbrt testCosmology testOutputList test27cellsStars \
+		 test27cellsStars_subset
 
 # Rebuild tests when SWIFT is updated.
 $(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -76,6 +78,12 @@ test27cells_subset_SOURCES = test27cells.c
 
 test27cells_subset_CFLAGS = $(AM_CFLAGS) -DTEST_DOSELF_SUBSET -DTEST_DOPAIR_SUBSET
 
+test27cellsStars_SOURCES = test27cellsStars.c
+
+test27cellsStars_subset_SOURCES = test27cellsStars.c
+
+test27cellsStars_subset_CFLAGS = $(AM_CFLAGS) -DTEST_DOSELF_SUBSET -DTEST_DOPAIR_SUBSET
+
 testPeriodicBC_SOURCES = testPeriodicBC.c
 
 test125cells_SOURCES = test125cells.c
@@ -130,4 +138,6 @@ EXTRA_DIST = testReading.sh makeInput.py testActivePair.sh \
 	     fft_params.yml tolerance_periodic_BC_normal.dat tolerance_periodic_BC_perturbed.dat \
 	     testEOS.sh testEOS_plot.sh testSelectOutput.sh selectOutput.yml \
              output_list_params.yml output_list_time.txt output_list_redshift.txt \
-             output_list_scale_factor.txt
+             output_list_scale_factor.txt testEOS.sh testEOS_plot.sh \
+	     test27cellsStars.sh test27cellsStarsPerturbed.sh star_tolerance_27_normal.dat \
+	     star_tolerance_27_perturbed.dat

tests/star_tolerance_27_normal.dat

+#   ID      pos_x      pos_y      pos_z     wcount     wcount_dh
+    0	    1e-6       1e-6	  1e-6 	    4e-4       1.2e-2
+    0	    1e-6       1e-6	  1e-6 	    1e-4       2e-4
+    0	    1e-6       1e-6	  1e-6 	    1e-6       1e-6

tests/star_tolerance_27_perturbed.dat

+#   ID      pos_x      pos_y      pos_z     wcount     wcount_dh
+    0	    1e-6       1e-6	  1e-6 	    2e-4       1e-2
+    0	    1e-6       1e-6	  1e-6 	    1e-5       2e-3
+    0	    1e-6       1e-6	  1e-6 	    1e-6       1e0

tests/test27cellsStars.c

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2015 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/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <fenv.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+
+/* Local headers. */
+#include "swift.h"
+
+#define DOSELF1 runner_doself_branch_star_density
+#define DOSELF1_SUBSET runner_doself_subset_branch_star_density
+#ifdef TEST_DOSELF_SUBSET
+#define DOSELF1_NAME "runner_doself_subset_branch_star_density"
+#else
+#define DOSELF1_NAME "runner_doself1_branch_star_density"
+#endif
+
+#define DOPAIR1 runner_dopair_branch_star_density
+#define DOPAIR1_SUBSET runner_dopair_subset_branch_star_density
+#ifdef TEST_DOPAIR_SUBSET
+#define DOPAIR1_NAME "runner_dopair_subset_branch_star_density"
+#else
+#define DOPAIR1_NAME "runner_dopair_branch_star_density"
+#endif
+
+#define NODE_ID 0
+
+/**
+ * @brief Constructs a cell and all of its particle in a valid state prior to
+ * a DOPAIR or DOSELF calcuation.
+ *
+ * @param n The cube root of the number of particles.
+ * @param n The cube root of the number of star particles.
+ * @param offset The position of the cell offset from (0,0,0).
+ * @param size The cell size.
+ * @param h The smoothing length of the particles in units of the inter-particle
+ * separation.
+ * @param partId The running counter of IDs.
+ * @param pert The perturbation to apply to the particles in the cell in units
+ * of the inter-particle separation.
+ * @param h_pert The perturbation to apply to the smoothing length.
+ */
+struct cell *make_cell(size_t n, size_t n_star, double *offset, double size, double h,
+                       long long *partId, long long *spartId, double pert, double h_pert) {
+  const size_t count = n * n * n;
+  const size_t scount = n_star * n_star * n_star;
+  float h_max = 0.f;
+  struct cell *cell = (struct cell *)malloc(sizeof(struct cell));
+  bzero(cell, sizeof(struct cell));
+
+  if (posix_memalign((void **)&cell->parts, part_align,
+                     count * sizeof(struct part)) != 0) {
+    error("couldn't allocate particles, no. of particles: %d", (int)count);
+  }
+  bzero(cell->parts, count * sizeof(struct part));
+
+  /* Construct the parts */
+  struct part *part = cell->parts;
+  for (size_t x = 0; x < n; ++x) {
+    for (size_t y = 0; y < n; ++y) {
+      for (size_t z = 0; z < n; ++z) {
+        part->x[0] =
+            offset[0] +
+            size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
+        part->x[1] =
+            offset[1] +
+            size * (y + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
+        part->x[2] =
+            offset[2] +
+            size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
+
+	part->v[0] = 0;
+	part->v[1] = 0;
+	part->v[2] = 0;
+	if (h_pert)
+	  part->h = size * h * random_uniform(1.f, h_pert) / (float)n;
+	else
+	  part->h = size * h / (float)n;
+	h_max = fmaxf(h_max, part->h);
+	part->id = ++(*partId);
+	
+	part->time_bin = 1;
+
+#ifdef SWIFT_DEBUG_CHECKS
+	part->ti_drift = 8;
+	part->ti_kick = 8;
+#endif
+	++part;
+      }
+    }
+  }
+
+  /* Construct the sparts */
+  if (posix_memalign((void **)&cell->sparts, part_align,
+                     scount * sizeof(struct spart)) != 0) {
+    error("couldn't allocate particles, no. of particles: %d", (int)scount);
+  }
+  bzero(cell->sparts, scount * sizeof(struct spart));
+
+  struct spart *spart = cell->sparts;
+  for (size_t x = 0; x < n_star; ++x) {
+    for (size_t y = 0; y < n_star; ++y) {
+      for (size_t z = 0; z < n_star; ++z) {
+        spart->x[0] =
+            offset[0] +
+            size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n_star;
+        spart->x[1] =
+            offset[1] +
+            size * (y + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n_star;
+        spart->x[2] =
+            offset[2] +
+            size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n_star;
+
+	spart->v[0] = 0;
+	spart->v[1] = 0;
+	spart->v[2] = 0;
+	if (h_pert)
+	  spart->h = size * h * random_uniform(1.f, h_pert) / (float)n_star;
+	else
+	  spart->h = size * h / (float)n_star;
+	h_max = fmaxf(h_max, spart->h);
+	spart->id = ++(*spartId);
+
+	spart->time_bin = 1;
+
+#ifdef SWIFT_DEBUG_CHECKS
+	spart->ti_drift = 8;
+	spart->ti_kick = 8;
+#endif
+
+	++spart;
+      }
+    }
+  }
+
+  /* Cell properties */
+  cell->split = 0;
+  cell->h_max = h_max;
+  cell->count = count;
+  cell->scount = scount;
+  cell->dx_max_part = 0.;
+  cell->dx_max_sort = 0.;
+  cell->width[0] = size;
+  cell->width[1] = size;
+  cell->width[2] = size;
+  cell->loc[0] = offset[0];
+  cell->loc[1] = offset[1];
+  cell->loc[2] = offset[2];
+
+  cell->ti_old_part = 8;
+  cell->ti_hydro_end_min = 8;
+  cell->ti_hydro_end_max = 8;
+  cell->nodeID = NODE_ID;
+
+  shuffle_particles(cell->parts, cell->count);
+  shuffle_sparticles(cell->sparts, cell->scount);
+
+  cell->sorted = 0;
+  for (int k = 0; k < 13; k++) cell->sort[k] = NULL;
+
+  return cell;
+}
+
+void clean_up(struct cell *ci) {
+  free(ci->parts);
+  free(ci->sparts);
+  for (int k = 0; k < 13; k++)
+    if (ci->sort[k] != NULL) free(ci->sort[k]);
+  free(ci);
+}
+
+/**
+ * @brief Initializes all particles field to be ready for a density calculation
+ */
+void zero_particle_fields(struct cell *c) {
+  for (int pid = 0; pid < c->scount; pid++) {
+    star_init_spart(&c->sparts[pid]);
+  }
+}
+
+/**
+ * @brief Ends the loop by adding the appropriate coefficients
+ */
+void end_calculation(struct cell *c, const struct cosmology *cosmo) {
+  for (int pid = 0; pid < c->scount; pid++) {
+    star_end_density(&c->sparts[pid], cosmo);
+
+    /* Recover the common "Neighbour number" definition */
+    c->sparts[pid].wcount *= pow_dimension(c->sparts[pid].h);
+    c->sparts[pid].wcount *= kernel_norm;
+  }
+}
+
+/**
+ * @brief Dump all the particles to a file
+ */
+void dump_particle_fields(char *fileName, struct cell *main_cell,
+                          struct cell **cells) {
+  FILE *file = fopen(fileName, "w");
+
+  /* Write header */
+  fprintf(file,
+          "# %4s %10s %10s %10s %13s %13s\n",
+          "ID", "pos_x", "pos_y", "pos_z", "wcount", "wcount_dh");
+
+  fprintf(file, "# Main cell --------------------------------------------\n");
+
+  /* Write main cell */
+  for (int pid = 0; pid < main_cell->scount; pid++) {
+    fprintf(file,
+            "%6llu %10f %10f %10f %13e %13e\n",
+            main_cell->sparts[pid].id, main_cell->sparts[pid].x[0],
+            main_cell->sparts[pid].x[1], main_cell->sparts[pid].x[2],
+            main_cell->sparts[pid].wcount,
+            main_cell->sparts[pid].wcount_dh);
+  }
+
+  /* Write all other cells */
+  for (int i = 0; i < 3; ++i) {
+    for (int j = 0; j < 3; ++j) {
+      for (int k = 0; k < 3; ++k) {
+        struct cell *cj = cells[i * 9 + j * 3 + k];
+        if (cj == main_cell) continue;
+
+        fprintf(file,
+                "# Offset: [%2d %2d %2d] -----------------------------------\n",
+                i - 1, j - 1, k - 1);
+
+        for (int pjd = 0; pjd < cj->scount; pjd++) {
+          fprintf(
+              file,
+              "%6llu %10f %10f %10f %13e %13e\n",
+              cj->sparts[pjd].id, cj->sparts[pjd].x[0], cj->sparts[pjd].x[1],
+              cj->sparts[pjd].x[2],
+              cj->sparts[pjd].wcount, cj->sparts[pjd].wcount_dh);
+        }