Pre-merge fixes

examples/EAGLE_12/eagle_12.yml

@@ -23,8 +23,9 @@ TimeIntegration:
   dt_max:     1e-4  # The maximal time-step size of the simulation (in internal units).
   
 Scheduler:
-  max_top_level_cells:    15
-  
+  max_top_level_cells:    8
+  cell_split_size:        64
+
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
@@ -42,9 +43,10 @@ Statistics:
 # Parameters for the self-gravity scheme
 Gravity:
   eta:                    0.025     # Constant dimensionless multiplier for time integration.
-  theta:                  0.85      # Opening angle (Multipole acceptance criterion)
+  theta:                  0.7       # Opening angle (Multipole acceptance criterion)
   comoving_softening:     0.0026994 # Comoving softening length (in internal units).
   max_physical_softening: 0.0007    # Physical softening length (in internal units).
+  mesh_side_length:       32
   
 # Parameters for the hydrodynamics scheme
 SPH:

examples/EAGLE_25/eagle_25.yml

@@ -23,7 +23,8 @@ TimeIntegration:
   dt_max:     1e-4  # The maximal time-step size of the simulation (in internal units).
 
 Scheduler:
-  max_top_level_cells:    20
+  max_top_level_cells:    16
+  cell_split_size:        64
 
 # Parameters governing the snapshots
 Snapshots:
@@ -41,9 +42,10 @@ Statistics:
 # Parameters for the self-gravity scheme
 Gravity:
   eta:                    0.025    # Constant dimensionless multiplier for time integration. 
-  theta:                  0.85     # Opening angle (Multipole acceptance criterion)
+  theta:                  0.7     # Opening angle (Multipole acceptance criterion)
   comoving_softening:     0.0026994 # Comoving softening length (in internal units).
   max_physical_softening: 0.0007    # Physical softening length (in internal units).
+  mesh_side_length:       32
 
 # Parameters for the hydrodynamics scheme
 SPH:

examples/ZeldovichPancake_3D/zeldovichPancake.yml

@@ -46,7 +46,7 @@ Cosmology:
 Gravity:
   mesh_side_length:   16
   eta: 0.025
-  theta: 0.7
+  theta: 0.85
   r_cut_max: 5.
   comoving_softening: 0.0001
   max_physical_softening: 0.0001

src/gravity_cache.h

@@ -136,7 +136,7 @@ static INLINE void gravity_cache_init(struct gravity_cache *c,
 }
 
 /**
- * @param Zero all the output fields (acceleration and potential) of a
+ * @brief Zero all the output fields (acceleration and potential) of a
  * #gravity_cache.
  *
  * @param c The #gravity_cache to zero.
@@ -328,12 +328,16 @@ gravity_cache_populate_no_mpole(const timebin_t max_active_bin,
  * the multi-pole.
  *
  * @param max_active_bin The largest active bin in the current time-step.
+ * @param periodic Are we using periodic BCs ?
+ * @param dim The size of the simulation volume along each dimension.
  * @param c The #gravity_cache to fill.
  * @param gparts The #gpart array to read from.
  * @param gcount The number of particles to read.
  * @param gcount_padded The number of particle to read padded to the next
  * multiple of the vector length.
  * @param cell The cell we play with (to get reasonable padding positions).
+ * @param CoM The position of the multipole.
+ * @param r_max2 The square of the multipole radius.
  * @param grav_props The global gravity properties.
  */
 __attribute__((always_inline)) INLINE static void

src/kernel_gravity.h

@@ -26,8 +26,6 @@
 #include "inline.h"
 #include "minmax.h"
 
-//#define GADGET2_SOFTENING_CORRECTION
-
 #ifdef GADGET2_SOFTENING_CORRECTION
 /*! Conversion factor between Plummer softening and internal softening */
 #define kernel_gravity_softening_plummer_equivalent 2.8

src/mesh_gravity.c

@@ -118,52 +118,6 @@ __attribute__((always_inline)) INLINE static void CIC_set(
   mesh[row_major_id_periodic(i + 1, j + 1, k + 1, N)] += value * dx * dy * dz;
 }
 
-/**
- * @brief Assigns a given multipole to a density mesh using the CIC method.
- *
- * @param m The #multipole.
- * @param rho The density mesh.
- * @param N the size of the mesh along one axis.
- * @param fac The width of a mesh cell.
- * @param dim The dimensions of the simulation box.
- */
-INLINE static void multipole_to_mesh_CIC(const struct gravity_tensors* m,
-                                         double* rho, int N, double fac,
-                                         const double dim[3]) {
-  error("aa");
-  /* Box wrap the multipole's position */
-  const double CoM_x = box_wrap(m->CoM[0], 0., dim[0]);
-  const double CoM_y = box_wrap(m->CoM[1], 0., dim[1]);
-  const double CoM_z = box_wrap(m->CoM[2], 0., dim[2]);
-
-  /* Workout the CIC coefficients */
-  int i = (int)(fac * CoM_x);
-  if (i >= N) i = N - 1;
-  const double dx = fac * CoM_x - i;
-  const double tx = 1. - dx;
-
-  int j = (int)(fac * CoM_y);
-  if (j >= N) j = N - 1;
-  const double dy = fac * CoM_y - j;
-  const double ty = 1. - dy;
-
-  int k = (int)(fac * CoM_z);
-  if (k >= N) k = N - 1;
-  const double dz = fac * CoM_z - k;
-  const double tz = 1. - dz;
-
-#ifdef SWIFT_DEBUG_CHECKS
-  if (i < 0 || i >= N) error("Invalid multipole position in x");
-  if (j < 0 || j >= N) error("Invalid multipole position in y");
-  if (k < 0 || k >= N) error("Invalid multipole position in z");
-#endif
-
-  const double mass = m->m_pole.M_000;
-
-  /* CIC ! */
-  CIC_set(rho, N, i, j, k, tx, ty, tz, dx, dy, dz, mass);
-}
-
 /**
  * @brief Assigns a given #gpart to a density mesh using the CIC method.
  *
@@ -245,9 +199,9 @@ void mesh_to_gparts_CIC(struct gpart* gp, const double* pot, int N, double fac,
   const double tz = 1. - dz;
 
 #ifdef SWIFT_DEBUG_CHECKS
-  if (i < 0 || i >= N) error("Invalid multipole position in x");
-  if (j < 0 || j >= N) error("Invalid multipole position in y");
-  if (k < 0 || k >= N) error("Invalid multipole position in z");
+  if (i < 0 || i >= N) error("Invalid gpart position in x");
+  if (j < 0 || j >= N) error("Invalid gpart position in y");
+  if (k < 0 || k >= N) error("Invalid gpart position in z");
 #endif
 
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
@@ -329,10 +283,8 @@ void pm_mesh_compute_potential(struct pm_mesh* mesh, const struct engine* e) {
   const struct space* s = e->s;
   const double r_s = mesh->r_s;
   const double box_size = s->dim[0];
-  const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
   const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
 
-  if (cdim[0] != cdim[1] || cdim[0] != cdim[2]) error("Non-square mesh");
   if (r_s <= 0.) error("Invalid value of a_smooth");
 
   /* Some useful constants */
@@ -340,20 +292,6 @@ void pm_mesh_compute_potential(struct pm_mesh* mesh, const struct engine* e) {
   const int N_half = N / 2;
   const double cell_fac = N / box_size;
 
-/* Recover the list of top-level multipoles */
-/* const int nr_cells = s->nr_cells; */
-/* struct gravity_tensors* restrict multipoles = s->multipoles_top; */
-
-#ifdef SWIFT_DEBUG_CHECKS
-/* const struct cell* cells = s->cells_top; */
-/* const integertime_t ti_current = e->ti_current; */
-
-/* /\* Make sure everything has been drifted to the current point *\/ */
-/* for (int i = 0; i < nr_cells; ++i) */
-/*   if (cells[i].ti_old_multipole != ti_current) */
-/*     error("Top-level multipole %d not drifted", i); */
-#endif
-
   /* Use the memory allocated for the potential to temporarily store rho */
   double* restrict rho = mesh->potential;
   if (rho == NULL) error("Error allocating memory for density mesh");
@@ -371,11 +309,9 @@ void pm_mesh_compute_potential(struct pm_mesh* mesh, const struct engine* e) {
   fftw_plan inverse_plan = fftw_plan_dft_c2r_3d(
       N, N, N, frho, rho, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
 
-  /* Do a CIC mesh assignment of the multipoles */
-  /* for (int i = 0; i < nr_cells; ++i) */
-  /*   multipole_to_mesh_CIC(&multipoles[i], rho, N, cell_fac, dim); */
   bzero(rho, N * N * N * sizeof(double));
 
+  /* Do a CIC mesh assignment of the gparts */
   for (size_t i = 0; i < e->s->nr_gparts; ++i)
     gpart_to_mesh_CIC(&e->s->gparts[i], rho, N, cell_fac, dim);
 

src/runner_doiact_grav.h

@@ -1376,7 +1376,7 @@ static INLINE void runner_dopair_recursive_grav(struct runner *r,
   } else if (!ci->split && !cj->split) {
 
     /* We have two leaves. Go P-P. */
-    runner_dopair_grav_pp(r, ci, cj, /*symmetric*/ 1, /*allow_mpoles*/ 0);
+    runner_dopair_grav_pp(r, ci, cj, /*symmetric*/ 1, /*allow_mpoles*/ 1);
 
   } else {
 
@@ -1584,10 +1584,8 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
 
       /* Call the PM interaction fucntion on the active sub-cells of ci */
 
-      runner_dopair_recursive_grav_pm(r, ci, cj);
-      // runner_dopair_grav_mm(r, ci, cj);
-      //runner_dopair_grav_pp(r, ci, cj, 0, 1);
-
+      //runner_dopair_recursive_grav_pm(r, ci, cj);
+      runner_dopair_grav_mm(r, ci, cj);
     } /* We are in charge of this pair */
   }   /* Loop over top-level cells */