No tree-walk, only P-P interactions.

examples/ZeldovichPancake_3D/makeIC.py

@@ -24,7 +24,7 @@ from numpy import *
 
 # Parameters
 gamma = 5./3.          # Gas adiabatic index
-numPart_1D = 32        # Number of particles
+numPart_1D = 64        # Number of particles
 lambda_i = 1.975e24    # h^-1 m (= 64 h^-1 Mpc)
 x_min = -0.5 * lambda_i
 x_max = 0.5 * lambda_i
@@ -35,9 +35,9 @@ z_c = 1.
 z_i = 100.
 fileName = "zeldovichPancake.hdf5"
 
-Mpc_in_m = 3.086e22
-Msol_in_kg = 1.988435e30
-Gyr_in_s = 3.154e16
+Mpc_in_m = 3.085678e22
+Msol_in_kg = 1.989e30
+Gyr_in_s = 3.085678e19
 mH_in_kg = 1.6737236e-27
 k_in_J_K = 1.38064852e-23
 

examples/ZeldovichPancake_3D/zeldovichPancake.yml

 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.989e43 # 10^10 M_su in grams
-  UnitLength_in_cgs:   3.085678e24   # Mpc in centimeters
+  UnitMass_in_cgs:     1.989e43  # 10^10 M_su in grams
+  UnitLength_in_cgs:   3.085678e24  # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5   # km/s in centimeters per second
   UnitCurrent_in_cgs:  1   # Amperes
   UnitTemp_in_cgs:     1   # Kelvin
@@ -42,6 +42,7 @@ Cosmology:
   a_end: 1.0
 
 Gravity:
+  mesh_side_length:   12
   eta: 0.025
   theta: 0.9
   comoving_softening: 0.0001

src/engine.c

@@ -6257,7 +6257,7 @@ void engine_recompute_displacement_constraint(struct engine *e) {
   /* Mesh forces smoothing scale */
   float a_smooth;
   if ((e->policy & engine_policy_self_gravity) && e->s->periodic == 1)
-    a_smooth = e->gravity_properties->a_smooth * e->s->dim[0] / e->s->cdim[0];
+    a_smooth = e->mesh->a_smooth;
   else
     a_smooth = FLT_MAX;
 

src/gravity/Default/gravity_iact.h

@@ -47,6 +47,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
     float r2, float h2, float h_inv, float h_inv3, float mass, float *f_ij,
     float *pot_ij) {
 
+  /* *f_ij = 0.f; */
+  /* *pot_ij = 0.f; */
+  /* return; */
+
   /* Get the inverse distance */
   const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
 
@@ -92,6 +96,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
     float r2, float h2, float h_inv, float h_inv3, float mass, float rlr_inv,
     float *f_ij, float *pot_ij) {
 
+  /* *f_ij = 0.f; */
+  /* *pot_ij = 0.f; */
+  /* return; */
+
   /* Get the inverse distance */
   const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
   const float r = r2 * r_inv;
@@ -147,6 +155,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
     float r_x, float r_y, float r_z, float r2, float h, float h_inv,
     const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) {
 
+  *f_x = 0.f;
+  *f_y = 0.f;
+  *f_z = 0.f;
+  *pot = 0.f;
+  return;
+
 #if SELF_GRAVITY_MULTIPOLE_ORDER < 3
   float f_ij;
   runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,

src/gravity_cache.h

@@ -185,7 +185,7 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
     const float dy = y[i] - CoM[1];
     const float dz = z[i] - CoM[2];
     const float r2 = dx * dx + dy * dy + dz * dz;
-    use_mpole[i] = gravity_M2P_accept(r_max2, theta_crit2, r2);
+    use_mpole[i] = 0 * gravity_M2P_accept(r_max2, theta_crit2, r2);
   }
 
 #ifdef SWIFT_DEBUG_CHECKS

src/kernel_long_gravity.h

@@ -30,6 +30,8 @@
 /* Standard headers */
 #include <math.h>
 
+#define GADGET2_LONG_RANGE_CORRECTION
+
 /**
  * @brief Computes the long-range correction term for the potential calculation
  * coming from FFT.

src/multipole.h

@@ -1523,8 +1523,6 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
                                const struct gravity_props *props, int periodic,
                                const double dim[3]) {
 
-  return;
-
   /* Recover some constants */
   const float eps = props->epsilon_cur;
   const float eps_inv = props->epsilon_cur_inv;
@@ -2263,8 +2261,6 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
   gp->num_interacted += lb->num_interacted;
 #endif
 
-  return;
-  
   /* Local accumulator */
   double a_grav[3] = {0., 0., 0.};
   double pot = 0.;

src/runner_doiact_fft.c

@@ -37,7 +37,7 @@
 #include "space.h"
 #include "timers.h"
 
-#ifdef HAVE_FFTW
+#ifdef HAVE_FFTW2
 
 /**
  * @brief Returns 1D index of a 3D NxNxN array using row-major style.
@@ -441,7 +441,7 @@ void print_carray(fftw_complex* array, int N) {
  */
 void runner_do_grav_fft(struct runner* r, int timer) {
 
-#ifdef HAVE_FFTW
+#ifdef HAVE_FFTW2
 
   const struct engine* e = r->e;
   const struct space* s = e->s;

src/runner_doiact_grav.h

@@ -29,6 +29,9 @@
 #include "space_getsid.h"
 #include "timers.h"
 
+static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
+                                         struct cell *cj, int symmetric);
+
 /**
  * @brief Recursively propagate the multipoles down the tree by applying the
  * L2L and L2P kernels.
@@ -57,6 +60,8 @@ static INLINE void runner_do_grav_down(struct runner *r, struct cell *c,
 
   if (c->split) { /* Node case */
 
+    message("aa");
+
     /* Add the field-tensor to all the 8 progenitors */
     for (int k = 0; k < 8; ++k) {
       struct cell *cp = c->progeny[k];
@@ -73,7 +78,7 @@ static INLINE void runner_do_grav_down(struct runner *r, struct cell *c,
         struct grav_tensor shifted_tensor;
 
         /* If the tensor received any contribution, push it down */
-        if (c->multipole->pot.interacted) {
+        if (1 || c->multipole->pot.interacted) {
 
           /* Shift the field tensor */
           gravity_L2L(&shifted_tensor, &c->multipole->pot, cp->multipole->CoM,
@@ -129,7 +134,7 @@ static INLINE void runner_do_grav_down(struct runner *r, struct cell *c,
  * @param ci The #cell with field tensor to interact.
  * @param cj The #cell with the multipole.
  */
-static INLINE void runner_dopair_grav_mm(const struct runner *r,
+static INLINE void runner_dopair_grav_mm(struct runner *r,
                                          struct cell *restrict ci,
                                          struct cell *restrict cj) {
 
@@ -144,6 +149,8 @@ static INLINE void runner_dopair_grav_mm(const struct runner *r,
 
   TIMER_TIC;
 
+  // error("aa");
+
   /* Anything to do here? */
   if (!cell_is_active_gravity(ci, e) || ci->nodeID != engine_rank) return;
 
@@ -169,8 +176,10 @@ static INLINE void runner_dopair_grav_mm(const struct runner *r,
         cj->ti_old_multipole, cj->nodeID, ci->nodeID, e->ti_current);
 
   /* Let's interact at this level */
-  gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
-              cj->multipole->CoM, props, periodic, dim);
+  if (0)
+    gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
+                cj->multipole->CoM, props, periodic, dim);
+  runner_dopair_grav_pp(r, ci, cj, 0);
 
   TIMER_TOC(timer_dopair_grav_mm);
 }
@@ -185,6 +194,8 @@ static INLINE void runner_dopair_grav_pp_full(const struct engine *e,
 
   TIMER_TIC;
 
+  error("aa");
+
   /* Loop over all particles in ci... */
   for (int pid = 0; pid < gcount_i; pid++) {
 
@@ -278,6 +289,8 @@ static INLINE void runner_dopair_grav_pp_truncated(
     int gcount_padded_j, struct gpart *restrict gparts_i,
     struct gpart *restrict gparts_j) {
 
+  const float dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
+
   TIMER_TIC;
 
   /* Loop over all particles in ci... */
@@ -324,9 +337,14 @@ static INLINE void runner_dopair_grav_pp_truncated(
       const float mass_j = cj_cache->m[pjd];
 
       /* Compute the pairwise (square) distance. */
-      const float dx = x_i - x_j;
-      const float dy = y_i - y_j;
-      const float dz = z_i - z_j;
+      float dx = x_i - x_j;
+      float dy = y_i - y_j;
+      float dz = z_i - z_j;
+
+      dx = nearestf(dx, dim[0]);
+      dy = nearestf(dy, dim[1]);
+      dz = nearestf(dz, dim[2]);
+
       const float r2 = dx * dx + dy * dy + dz * dz;
 
 #ifdef SWIFT_DEBUG_CHECKS
@@ -449,7 +467,7 @@ static INLINE void runner_dopair_grav_pm(
  * @param cj The other #cell.
  */
 static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
-                                         struct cell *cj) {
+                                         struct cell *cj, int symmetric) {
 
   const struct engine *e = r->e;
 
@@ -468,10 +486,8 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
   /* Recover some useful constants */
   struct space *s = e->s;
   const int periodic = s->periodic;
-  const double cell_width = s->width[0];
-  const double a_smooth = e->gravity_properties->a_smooth;
   const double r_cut_min = e->gravity_properties->r_cut_min;
-  const double rlr = cell_width * a_smooth;
+  const double rlr = e->mesh->a_smooth;
   const double min_trunc = rlr * r_cut_min;
   const float rlr_inv = 1. / rlr;
 
@@ -481,7 +497,8 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
 
   /* Get the distance vector between the pairs, wrapping. */
   double cell_shift[3];
-  space_getsid(s, &ci, &cj, cell_shift);
+  struct cell *cii = ci, *cjj = cj;
+  space_getsid(s, &cii, &cjj, cell_shift);
 
   /* Record activity status */
   const int ci_active = cell_is_active_gravity(ci, e);
@@ -494,14 +511,15 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
     error("Un-drifted multipole");
 
   /* Centre of the cell pair */
-  const double loc[3] = {ci->loc[0],   // + 0. * ci->width[0],
-                         ci->loc[1],   // + 0. * ci->width[1],
-                         ci->loc[2]};  // + 0. * ci->width[2]};
+  /* const double loc[3] = {ci->loc[0],   // + 0. * ci->width[0], */
+  /*                        ci->loc[1],   // + 0. * ci->width[1], */
+  /*                        ci->loc[2]};  // + 0. * ci->width[2]}; */
 
   /* Shift to apply to the particles in each cell */
-  const double shift_i[3] = {loc[0] + cell_shift[0], loc[1] + cell_shift[1],
-                             loc[2] + cell_shift[2]};
-  const double shift_j[3] = {loc[0], loc[1], loc[2]};
+  const double shift_i[3] = {
+      0., 0., 0.};  //{loc[0] + cell_shift[0], loc[1] + cell_shift[1],
+                    // loc[2] + cell_shift[2]};
+  const double shift_j[3] = {0., 0., 0.};  //{loc[0], loc[1], loc[2]};
 
   /* Recover the multipole info and shift the CoM locations */
   const float rmax_i = ci->multipole->r_max;
@@ -556,7 +574,7 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
       runner_dopair_grav_pm(e, ci_cache, gcount_i, gcount_padded_i, ci->gparts,
                             CoM_j, multi_j, cj);
     }
-    if (cj_active) {
+    if (cj_active && symmetric) {
 
       /* First the P2P */
       runner_dopair_grav_pp_full(e, cj_cache, ci_cache, gcount_j, gcount_i,
@@ -652,6 +670,8 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
 static INLINE void runner_doself_grav_pp_full(struct runner *r,
                                               struct cell *c) {
 
+  error("bb");
+
   /* Some constants */
   const struct engine *const e = r->e;
   struct gravity_cache *const ci_cache = &r->ci_gravity_cache;
@@ -778,10 +798,7 @@ static INLINE void runner_doself_grav_pp_truncated(struct runner *r,
 
   /* Some constants */
   const struct engine *const e = r->e;
-  const struct space *s = e->s;
-  const double cell_width = s->width[0];
-  const double a_smooth = e->gravity_properties->a_smooth;
-  const double rlr = cell_width * a_smooth;
+  const double rlr = e->mesh->a_smooth;
   const float rlr_inv = 1. / rlr;
 
   /* Caches to play with */
@@ -908,10 +925,9 @@ static INLINE void runner_doself_grav_pp(struct runner *r, struct cell *c) {
   const struct engine *e = r->e;
   const struct space *s = e->s;
   const int periodic = s->periodic;
-  const double cell_width = s->width[0];
-  const double a_smooth = e->gravity_properties->a_smooth;
+  const double a_smooth = e->mesh->a_smooth;
   const double r_cut_min = e->gravity_properties->r_cut_min;
-  const double min_trunc = cell_width * r_cut_min * a_smooth;
+  const double min_trunc = r_cut_min * a_smooth;
 
   TIMER_TIC;
 
@@ -965,11 +981,10 @@ static INLINE void runner_dopair_grav(struct runner *r, struct cell *ci,
   const struct space *s = e->s;
   const int nodeID = e->nodeID;
   const int periodic = s->periodic;
-  const double cell_width = s->width[0];
   const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
   const struct gravity_props *props = e->gravity_properties;
   const double theta_crit2 = props->theta_crit2;
-  const double max_distance = props->a_smooth * props->r_cut_max * cell_width;
+  const double max_distance = e->mesh->a_smooth * props->r_cut_max;
   const double max_distance2 = max_distance * max_distance;
 
   /* Anything to do here? */
@@ -1039,7 +1054,7 @@ static INLINE void runner_dopair_grav(struct runner *r, struct cell *ci,
   }
   /* We have two leaves. Go P-P. */
   else if (!ci->split && !cj->split) {
-    runner_dopair_grav_pp(r, ci, cj);
+    runner_dopair_grav_pp(r, ci, cj, 1);
   }
   /* Alright, we'll have to split and recurse. */
   else {
@@ -1169,10 +1184,9 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
   const struct space *s = e->s;
   const struct gravity_props *props = e->gravity_properties;
   const int periodic = s->periodic;
-  const double cell_width = s->width[0];
   const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
   const double theta_crit2 = props->theta_crit2;
-  const double max_distance = props->a_smooth * props->r_cut_max * cell_width;
+  const double max_distance = e->mesh->a_smooth * props->r_cut_max;
   const double max_distance2 = max_distance * max_distance;
 
   TIMER_TIC;