Fix inversion arguments in star_density

configure.ac

@@ -1038,6 +1038,25 @@ case "$with_hydro" in
    ;;
 esac
 
+# Check if debugging interactions stars is switched on.
+AC_ARG_ENABLE([debug-interactions-stars],
+   [AS_HELP_STRING([--enable-debug-interactions-stars],
+     [Activate interaction debugging for stars, logging a maximum of @<:@N@:>@ neighbours. Defaults to 256 if no value set.]
+   )],
+   [enable_debug_interactions_stars="$enableval"],
+   [enable_debug_interactions_stars="no"]
+)
+if test "$enable_debug_interactions_stars" != "no"; then
+    AC_DEFINE([DEBUG_INTERACTIONS_STARS],1,[Enable interaction debugging for stars])
+    if test "$enable_debug_interactions_stars" == "yes"; then
+      AC_DEFINE([MAX_NUM_OF_NEIGHBOURS_STARS],256,[The maximum number of particle neighbours to be logged for stars])
+      [enable_debug_interactions_stars="yes (Logging up to 256 neighbours)"]
+    else
+      AC_DEFINE_UNQUOTED([MAX_NUM_OF_NEIGHBOURS_STARS], [$enableval] ,[The maximum number of particle neighbours to be logged for stars])
+      [enable_debug_interactions_stars="yes (Logging up to $enableval neighbours)"]
+    fi
+fi
+
 # Check if debugging interactions is switched on.
 AC_ARG_ENABLE([debug-interactions],
    [AS_HELP_STRING([--enable-debug-interactions],
@@ -1061,6 +1080,7 @@ if test "$enable_debug_interactions" != "no"; then
   fi
 fi
 
+
 # SPH Kernel function
 AC_ARG_WITH([kernel],
    [AS_HELP_STRING([--with-kernel=<kernel>],

src/runner_doiact_star.h

@@ -156,7 +156,7 @@ void runner_dosubpair_star_density(struct runner *r, struct cell *restrict ci,
 #endif
 
       if (r2 < hig2)
-	runner_iact_nonsym_star_density(r2, dx, hj, hi, si, pj, a, H);
+	runner_iact_nonsym_star_density(r2, dx, hi, hj, si, pj, a, H);
 
     } /* loop over the parts in cj. */
   }   /* loop over the parts in ci. */

src/stars/Default/star.h

@@ -55,8 +55,8 @@ __attribute__((always_inline)) INLINE static void star_first_init_spart(
 __attribute__((always_inline)) INLINE static void star_init_spart(
   struct spart* sp) {
 
-#ifdef DEBUG_INTERACTIONS_SPH
-  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS; ++i) sp->ids_ngbs_density[i] = -1;
+#ifdef DEBUG_INTERACTIONS_STARS
+  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i) sp->ids_ngbs_density[i] = -1;
   sp->num_ngb_density = 0;
 #endif
 
@@ -151,13 +151,6 @@ __attribute__((always_inline)) INLINE static void star_prepare_force(
  * @param sp The particle to act upon
  */
 __attribute__((always_inline)) INLINE static void star_reset_acceleration(
-    struct spart *restrict p) {
-
-#ifdef DEBUG_INTERACTIONS_SPH
-  for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS; ++i) sp->ids_ngbs_force[i] = -1;
-  sp->num_ngb_force = 0;
-#endif
-
-}
+    struct spart *restrict p) {}
 
 #endif /* SWIFT_DEFAULT_STAR_H */

src/stars/Default/star_iact.h

@@ -27,11 +27,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_star_densit
 
   /* Compute contribution to the number of neighbours */
   si->wcount += wi;
-  si->wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+  si->wcount_dh -= (hydro_dimension * wi + ui * wi_dx);    
 
 #ifdef DEBUG_INTERACTIONS_STARS
   /* Update ngb counters */
-  if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS)
+  if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_STARS)
     si->ids_ngbs_density[si->num_ngb_density] = pj->id;
   ++si->num_ngb_density;
 #endif

src/stars/Default/star_part.h

@@ -68,7 +68,7 @@ struct spart {
 
 #ifdef DEBUG_INTERACTIONS_STARS
   /*! List of interacting particles in the density SELF and PAIR */
-  long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS];
+  long long ids_ngbs_density[MAX_NUM_OF_NEIGHBOURS_STARS];
 
   /*! Number of interactions in the density SELF and PAIR */
   int num_ngb_density;

src/tools.c

@@ -369,7 +369,6 @@ void pairs_all_star_density(struct runner *r, struct cell *ci, struct cell *cj)
 
       /* Hit or miss? */
       if (r2 < hig2) {
-
         /* Interact */
         runner_iact_nonsym_star_density(r2, dx, hi, pj->h, spi, pj, a, H);
       }
@@ -400,7 +399,6 @@ void pairs_all_star_density(struct runner *r, struct cell *ci, struct cell *cj)
 
       /* Hit or miss? */
       if (r2 < hjg2) {
-
         /* Interact */
         runner_iact_nonsym_star_density(r2, dx, hj, pi->h, spj, pi, a, H);
       }
@@ -502,7 +500,7 @@ 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];
+  float r2, hi, hj, hig2, dxi[3];
   struct spart *spi;
   struct part *pj;
   const struct engine *e = r->e;
@@ -534,7 +532,6 @@ void self_all_star_density(struct runner *r, struct cell *ci) {
 
       /* Hit or miss? */
       if (r2 > 0.f && r2 < hig2) {
-
         /* Interact */
         runner_iact_nonsym_star_density(r2, dxi, hi, hj, spi, pj, a, H);
       }

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-4       2e-3
     0	    1e-6       1e-6	  1e-6 	    1e-6       1e-6

tests/star_tolerance_27_perturbed_h.dat

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

tests/star_tolerance_27_perturbed_h2.dat

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

tests/test27cellsStars.c

@@ -114,7 +114,7 @@ struct cell *make_cell(size_t n, size_t n_star, double *offset, double size, dou
   }
 
   /* Construct the sparts */
-  if (posix_memalign((void **)&cell->sparts, part_align,
+  if (posix_memalign((void **)&cell->sparts, spart_align,
                      scount * sizeof(struct spart)) != 0) {
     error("couldn't allocate particles, no. of particles: %d", (int)scount);
   }
@@ -150,7 +150,6 @@ struct cell *make_cell(size_t n, size_t n_star, double *offset, double size, dou
 	spart->ti_drift = 8;
 	spart->ti_kick = 8;
 #endif
-
 	++spart;
       }
     }
@@ -401,13 +400,14 @@ int main(int argc, char *argv[]) {
   struct cell *cells[27];
   struct cell *main_cell;
   static long long partId = 0;
-  long long spartId = particles * particles * particles * 27 ;
+  long long spartId = particles * particles * particles * 27;
+
   for (int i = 0; i < 3; ++i) {
     for (int j = 0; j < 3; ++j) {
       for (int k = 0; k < 3; ++k) {
         double offset[3] = {i * size, j * size, k * size};
         cells[i * 9 + j * 3 + k] =
-	  make_cell(particles, sparticles, offset, size, h, &partId, &spartId + 1, perturbation, h_pert);
+	  make_cell(particles, sparticles, offset, size, h, &partId, &spartId, perturbation, h_pert);
 
         runner_do_drift_part(&runner, cells[i * 9 + j * 3 + k], 0);