diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 66a475f32ec06eb40ff2bc890bc156f76e3b7b9f..ad927ed2c55e008af2a98529bf1bd2d175a75f33 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -175,6 +175,11 @@ __attribute__((always_inline)) INLINE static void hydro_timestep_extra(
__attribute__((always_inline)) INLINE static void hydro_init_part(
struct part *restrict p, const struct hydro_space *hs) {
+#ifdef DEBUG_INTERACTIONS
+ for (int i = 0; i < NUM_OF_NEIGHBOURS; ++i) p->ids_ngbs_density[i] = -1;
+ p->num_ngb_density = 0;
+#endif
+
p->rho = 0.f;
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
@@ -311,6 +316,11 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
struct part *restrict p) {
+#ifdef DEBUG_INTERACTIONS
+ for (int i = 0; i < NUM_OF_NEIGHBOURS; ++i) p->ids_ngbs_force[i] = -1;
+ p->num_ngb_force = 0;
+#endif
+
/* Reset the acceleration. */
p->a_hydro[0] = 0.0f;
p->a_hydro[1] = 0.0f;
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 45ee89ceb1ba119b396aa847b39a8b6ff1b38b8c..a8affdc6e13b4d7e31416a365d6a6d0e7443d42c 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -103,6 +103,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pj->density.rot_v[0] += facj * curlvr[0];
pj->density.rot_v[1] += facj * curlvr[1];
pj->density.rot_v[2] += facj * curlvr[2];
+
+ /* Update ngb counters */
+#ifdef DEBUG_INTERACTIONS
+ if(pi->id == CHECK_PART_ID) pi->ids_ngbs_density[pi->num_ngb_density] = pj->id;
+ ++pi->num_ngb_density;
+ if(pj->id == CHECK_PART_ID) pj->ids_ngbs_density[pj->num_ngb_density] = pi->id;
+ ++pj->num_ngb_density;
+#endif
+
}
/**
@@ -151,6 +160,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
pi->density.rot_v[0] += fac * curlvr[0];
pi->density.rot_v[1] += fac * curlvr[1];
pi->density.rot_v[2] += fac * curlvr[2];
+
+#ifdef DEBUG_INTERACTIONS
+ /* Update ngb counters */
+ if(pi->id == CHECK_PART_ID) pi->ids_ngbs_density[pi->num_ngb_density] = pj->id;
+ ++pi->num_ngb_density;
+#endif
+
}
#ifdef WITH_VECTORIZATION
@@ -480,6 +496,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Change in entropy */
pi->entropy_dt += mj * visc_term * dvdr;
pj->entropy_dt += mi * visc_term * dvdr;
+
+#ifdef DEBUG_INTERACTIONS
+ /* Update ngb counters */
+ if(pi->id == CHECK_PART_ID) {
+ pi->ids_ngbs_force[pi->num_ngb_force] = pj->id;
+ }
+ ++pi->num_ngb_force;
+ if(pj->id == CHECK_PART_ID) {
+ pj->ids_ngbs_force[pj->num_ngb_force] = pi->id;
+ }
+ ++pj->num_ngb_force;
+#endif
+
}
/**
@@ -569,6 +598,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Change in entropy */
pi->entropy_dt += mj * visc_term * dvdr;
+
+#ifdef DEBUG_INTERACTIONS
+ /* Update ngb counters */
+ if(pi->id == CHECK_PART_ID) {
+ pi->ids_ngbs_force[pi->num_ngb_force] = pj->id;
+ }
+ ++pi->num_ngb_force;
+#endif
+
}
#ifdef WITH_VECTORIZATION
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 3e46b2351eb6a3871946dd9c69c8d108b10da0df..f4b1cca41d80a8efcc6c9a53adfb6d683383b8e3 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -77,6 +77,10 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
*num_fields = 10;
+#ifdef DEBUG_INTERACTIONS
+ *num_fields = 14;
+#endif
+
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
parts, x);
@@ -99,6 +103,17 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
parts, rho, convert_u);
list[9] = io_make_output_field_convert_part(
"Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, rho, convert_P);
+#ifdef DEBUG_INTERACTIONS
+ list[10] = io_make_output_field("Num_ngb_density", INT, 1, UNIT_CONV_NO_UNITS,
+ parts, num_ngb_density);
+ list[11] = io_make_output_field("Num_ngb_force", INT, 1, UNIT_CONV_NO_UNITS,
+ parts, num_ngb_force);
+ list[12] = io_make_output_field("Ids_ngb_density", LONGLONG, NUM_OF_NEIGHBOURS,
+ UNIT_CONV_NO_UNITS, parts, ids_ngbs_density);
+ list[13] = io_make_output_field("Ids_ngb_force", LONGLONG, NUM_OF_NEIGHBOURS,
+ UNIT_CONV_NO_UNITS, parts, ids_ngbs_force);
+#endif
+
}
/**
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 571aaf39ed2509d95e9f68c34ab8e6c09ded3fbb..52ff778f7e091e0acc5b215517378f2480409f55 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -33,6 +33,10 @@
#include "cooling_struct.h"
+#define NUM_OF_NEIGHBOURS 128
+#define CHECK_PART_ID 999999999999999999999
+//#define CHECK_PART_ID 5673486482283
+
/* Extra particle data not needed during the SPH loops over neighbours. */
struct xpart {
@@ -143,6 +147,20 @@ struct part {
#endif
+#ifdef DEBUG_INTERACTIONS
+ /*! List of interacting particles in the density SELF and PAIR */
+ long long ids_ngbs_density[NUM_OF_NEIGHBOURS];
+
+ /*! List of interacting particles in the force SELF and PAIR */
+ long long ids_ngbs_force[NUM_OF_NEIGHBOURS];
+
+ /*! Number of interactions in the density SELF and PAIR */
+ int num_ngb_density;
+
+ /*! Number of interactions in the force SELF and PAIR */
+ int num_ngb_force;
+#endif
+
} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_GADGET2_HYDRO_PART_H */
diff --git a/src/runner_doiact_vec.c b/src/runner_doiact_vec.c
index 8188663119351d88abdf3593f5fcb92a21645295..eecfa75c34bda78ba7a79617974f43b042b7c310 100644
--- a/src/runner_doiact_vec.c
+++ b/src/runner_doiact_vec.c
@@ -662,6 +662,26 @@ __attribute__((always_inline)) INLINE void runner_doself1_density_vec(
doi_mask = doi_mask & doi_mask_self_check;
doi_mask2 = doi_mask2 & doi_mask2_self_check;
+#ifdef DEBUG_INTERACTIONS
+ for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+ if (doi_mask & (1 << bit_index)) {
+ if(pi->id == CHECK_PART_ID) {
+ pi->ids_ngbs_density[pi->num_ngb_density] = parts[pjd + bit_index].id;
+ }
+ ++pi->num_ngb_density;
+ }
+ }
+
+ for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+ if (doi_mask2 & (1 << bit_index)) {
+ if(pi->id == CHECK_PART_ID) {
+ pi->ids_ngbs_density[pi->num_ngb_density] = parts[pjd + VEC_SIZE + bit_index].id;
+ }
+ ++pi->num_ngb_density;
+ }
+ }
+#endif
+
/* If there are any interactions left pack interaction values into c2
* cache. */
if (doi_mask) {
@@ -776,6 +796,12 @@ __attribute__((always_inline)) INLINE void runner_doself2_force_vec(
/* Get a pointer to the ith particle. */
pi = &parts[pid];
+
+ if(pi->id == CHECK_PART_ID) {
+ message("Here");
+ }
+
+
/* Is the ith particle active? */
if (!part_is_active_no_debug(pi, max_active_bin)) continue;
@@ -881,6 +907,17 @@ __attribute__((always_inline)) INLINE void runner_doself2_force_vec(
vec_combine_masks(v_doi_mask, v_doi_mask_self_check);
doi_mask = vec_form_int_mask(v_doi_mask);
+#ifdef DEBUG_INTERACTIONS
+ for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+ if (doi_mask & (1 << bit_index)) {
+ if(pi->id == CHECK_PART_ID) {
+ pi->ids_ngbs_force[pi->num_ngb_force] = parts[pjd + bit_index].id;
+ }
+ ++pi->num_ngb_force;
+ }
+ }
+#endif
+
/* If there are any interactions perform them. */
if (doi_mask) {
vector v_hj_inv;
@@ -1171,6 +1208,17 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
/* Form integer mask. */
doi_mask = vec_form_int_mask(v_doi_mask);
+#ifdef DEBUG_INTERACTIONS
+ for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+ if (doi_mask & (1 << bit_index)) {
+ if(pi->id == CHECK_PART_ID) {
+ pi->ids_ngbs_density[pi->num_ngb_density] = parts_j[sort_j[pjd + bit_index].i].id;
+ }
+ ++pi->num_ngb_density;
+ }
+ }
+#endif
+
/* If there are any interactions perform them. */
if (doi_mask)
runner_iact_nonsym_1_vec_density(
@@ -1302,6 +1350,17 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
/* Form integer mask. */
doj_mask = vec_form_int_mask(v_doj_mask);
+#ifdef DEBUG_INTERACTIONS
+ for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+ if (doj_mask & (1 << bit_index)) {
+ if(pj->id == CHECK_PART_ID) {
+ pj->ids_ngbs_density[pj->num_ngb_density] = parts_i[sort_i[ci_cache_idx + first_pi_align + bit_index].i].id;
+ }
+ ++pj->num_ngb_density;
+ }
+ }
+#endif
+
/* If there are any interactions perform them. */
if (doj_mask)
runner_iact_nonsym_1_vec_density(
@@ -1509,6 +1568,10 @@ void runner_dopair2_force_vec(struct runner *r, struct cell *ci,
struct part *restrict pi = &parts_i[sort_i[pid].i];
if (!part_is_active(pi, e)) continue;
+ if(pi->id == CHECK_PART_ID) {
+ message("Here");
+ }
+
/* Set the cache index. */
int ci_cache_idx = pid - first_pi_align;
@@ -1612,6 +1675,17 @@ void runner_dopair2_force_vec(struct runner *r, struct cell *ci,
/* Form integer masks. */
doi_mask = vec_form_int_mask(v_doi_mask);
+#ifdef DEBUG_INTERACTIONS
+ for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+ if (doi_mask & (1 << bit_index)) {
+ if(pi->id == CHECK_PART_ID) {
+ pi->ids_ngbs_force[pi->num_ngb_force] = parts_j[sort_j[pjd + bit_index].i].id;
+ }
+ ++pi->num_ngb_force;
+ }
+ }
+#endif
+
/* If there are any interactions perform them. */
if (doi_mask) {
vector v_hj_inv;
@@ -1652,6 +1726,10 @@ void runner_dopair2_force_vec(struct runner *r, struct cell *ci,
struct part *restrict pj = &parts_j[sort_j[pjd].i];
if (!part_is_active(pj, e)) continue;
+ if(pj->id == CHECK_PART_ID) {
+ message("Here");
+ }
+
/* Set the cache index. */
int cj_cache_idx = pjd;
@@ -1659,7 +1737,7 @@ void runner_dopair2_force_vec(struct runner *r, struct cell *ci,
/* Skip this particle if no particle in ci is within range of it. */
const float hj = cj_cache->h[cj_cache_idx];
const double dj_test =
- sort_j[pjd].d - hj * kernel_gamma - dx_max - rshift;
+ sort_j[pjd].d - hj * kernel_gamma - dx_max;
if (dj_test > di_max) continue;
/* Determine the exit iteration of the interaction loop. */
@@ -1755,6 +1833,17 @@ void runner_dopair2_force_vec(struct runner *r, struct cell *ci,
/* Form integer masks. */
doj_mask = vec_form_int_mask(v_doj_mask);
+#ifdef DEBUG_INTERACTIONS
+ for (int bit_index = 0; bit_index < VEC_SIZE; bit_index++) {
+ if (doj_mask & (1 << bit_index)) {
+ if(pj->id == CHECK_PART_ID) {
+ pj->ids_ngbs_force[pj->num_ngb_force] = parts_i[sort_i[ci_cache_idx + first_pi_align + bit_index].i].id;
+ }
+ ++pj->num_ngb_force;
+ }
+ }
+#endif
+
/* If there are any interactions perform them. */
if (doj_mask) {
vector v_hi_inv;