Do not access the space any more in the gravity calculation functions. All...

Do not access the space any more in the gravity calculation functions. All relevant information are stored in the pm_mesh structure.

examples/main.c

@@ -730,12 +730,14 @@ int main(int argc, char *argv[]) {
     /* Initialise the long-range gravity mesh */
     if (with_self_gravity && periodic) {
 #ifdef HAVE_FFTW
-      pm_mesh_init(&mesh, &gravity_properties, dim[0]);
+      pm_mesh_init(&mesh, &gravity_properties, dim);
 #else
       /* Need the FFTW library if periodic and self gravity. */
       error(
           "No FFTW library found. Cannot compute periodic long-range forces.");
 #endif
+    } else {
+      pm_mesh_init_no_mesh(&mesh, dim);
     }
 
     /* Initialize the space with these data. */

src/gravity_properties.c

@@ -53,10 +53,10 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
   p->mesh_size = parser_get_param_int(params, "Gravity:mesh_side_length");
   p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth",
                                            gravity_props_default_a_smooth);
-  p->r_cut_max = parser_get_opt_param_float(params, "Gravity:r_cut_max",
-                                            gravity_props_default_r_cut_max);
-  p->r_cut_min = parser_get_opt_param_float(params, "Gravity:r_cut_min",
-                                            gravity_props_default_r_cut_min);
+  p->r_cut_max_ratio = parser_get_opt_param_float(
+      params, "Gravity:r_cut_max", gravity_props_default_r_cut_max);
+  p->r_cut_min_ratio = parser_get_opt_param_float(
+      params, "Gravity:r_cut_min", gravity_props_default_r_cut_min);
 
   if (p->mesh_size % 2 != 0)
     error("The mesh side-length must be an even number.");
@@ -130,8 +130,9 @@ void gravity_props_print(const struct gravity_props *p) {
   message("Self-gravity mesh side-length: N=%d", p->mesh_size);
   message("Self-gravity mesh smoothing-scale: a_smooth=%f", p->a_smooth);
 
-  message("Self-gravity tree cut-off ratio: r_cut_max=%f", p->r_cut_max);
-  message("Self-gravity truncation cut-off ratio: r_cut_min=%f", p->r_cut_min);
+  message("Self-gravity tree cut-off ratio: r_cut_max=%f", p->r_cut_max_ratio);
+  message("Self-gravity truncation cut-off ratio: r_cut_min=%f",
+          p->r_cut_min_ratio);
 
   message("Self-gravity mesh truncation function: %s",
           kernel_long_gravity_truncation_name);
@@ -161,8 +162,8 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
   io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
   io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
   io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
-  io_write_attribute_f(h_grpgrav, "Mesh r_cut_max ratio", p->r_cut_max);
-  io_write_attribute_f(h_grpgrav, "Mesh r_cut_min ratio", p->r_cut_min);
+  io_write_attribute_f(h_grpgrav, "Mesh r_cut_max ratio", p->r_cut_max_ratio);
+  io_write_attribute_f(h_grpgrav, "Mesh r_cut_min ratio", p->r_cut_min_ratio);
   io_write_attribute_f(h_grpgrav, "Tree update frequency",
                        p->rebuild_frequency);
   io_write_attribute_s(h_grpgrav, "Mesh truncation function",

src/gravity_properties.h

@@ -47,11 +47,11 @@ struct gravity_props {
 
   /*! Distance below which the truncated mesh force is Newtonian in units of
    * a_smooth */
-  float r_cut_min;
+  float r_cut_min_ratio;
 
   /*! Distance above which the truncated mesh force is negligible in units of
    * a_smooth */
-  float r_cut_max;
+  float r_cut_max_ratio;
 
   /*! Time integration dimensionless multiplier */
   float eta;

src/mesh_gravity.c

@@ -508,19 +508,29 @@ void pm_mesh_interpolate_forces(const struct pm_mesh* mesh,
  *
  * @param mesh The #pm_mesh to initialise.
  * @param props The propoerties of the gravity scheme.
- * @param box_size The (comoving) side-length of the simulation volume.
+ * @param dim The (comoving) side-lengths of the simulation volume.
  */
 void pm_mesh_init(struct pm_mesh* mesh, const struct gravity_props* props,
-                  double box_size) {
+                  double dim[3]) {
 
 #ifdef HAVE_FFTW
 
+  if (dim[0] != dim[1] || dim[0] != dim[1])
+    error("Doing mesh-gravity on a non-cubic domain");
+
   const int N = props->mesh_size;
+  const double box_size = dim[0];
+
+  mesh->periodic = 1;
   mesh->N = N;
-  mesh->box_size = box_size;
+  mesh->dim[0] = dim[0];
+  mesh->dim[1] = dim[1];
+  mesh->dim[2] = dim[2];
   mesh->cell_fac = N / box_size;
   mesh->r_s = props->a_smooth * box_size / N;
   mesh->r_s_inv = 1. / mesh->r_s;
+  mesh->r_cut_max = mesh->r_s * props->r_cut_max_ratio;
+  mesh->r_cut_min = mesh->r_s * props->r_cut_min_ratio;
 
   /* Allocate the memory for the combined density and potential array */
   mesh->potential = (double*)fftw_malloc(sizeof(double) * N * N * N);
@@ -531,6 +541,30 @@ void pm_mesh_init(struct pm_mesh* mesh, const struct gravity_props* props,
 #endif
 }
 
+/**
+ * @brief Initialises the mesh for the case where we don't do mesh gravity
+ * calculations
+ *
+ * Crucially this set the 'periodic' propoerty to 0 and all the relevant values
+ * to a
+ * state where all calculations will default to pure non-periodic Newtonian.
+ *
+ * @param mesh The #pm_mesh to initialise.
+ * @param dim The (comoving) side-lengths of the simulation volume.
+ */
+void pm_mesh_init_no_mesh(struct pm_mesh* mesh, double dim[3]) {
+
+  bzero(mesh, sizeof(struct pm_mesh));
+
+  /* Fill in non-zero properties */
+  mesh->dim[0] = dim[0];
+  mesh->dim[1] = dim[1];
+  mesh->dim[2] = dim[2];
+  mesh->r_s = FLT_MAX;
+  mesh->r_cut_min = FLT_MAX;
+  mesh->r_cut_max = FLT_MAX;
+}
+
 /**
  * @brief Frees the memory allocated for the long-range mesh.
  */

src/mesh_gravity.h

@@ -34,14 +34,17 @@ struct gpart;
  */
 struct pm_mesh {
 
+  /*! Is the calculation using periodic BCs? */
+  int periodic;
+
   /*! Side-length of the mesh */
   int N;
 
   /*! Conversion factor between box and mesh size */
   double cell_fac;
 
-  /*! (Comoving) side-length of the volume covered by the mesh */
-  double box_size;
+  /*! (Comoving) side-length of the box along the three axis */
+  double dim[3];
 
   /*! Scale over which we smooth the forces */
   double r_s;
@@ -49,12 +52,19 @@ struct pm_mesh {
   /*! Inverse of the scale over which we smooth the forces */
   double r_s_inv;
 
+  /*! Distance beyond which tree forces are neglected */
+  double r_cut_max;
+
+  /*! Distance below which tree forces are Newtonian */
+  double r_cut_min;
+
   /*! Potential field */
   double *potential;
 };
 
 void pm_mesh_init(struct pm_mesh *mesh, const struct gravity_props *props,
-                  double box_size);
+                  double dim[3]);
+void pm_mesh_init_no_mesh(struct pm_mesh *mesh, double dim[3]);
 void pm_mesh_compute_potential(struct pm_mesh *mesh, const struct engine *e);
 void pm_mesh_interpolate_forces(const struct pm_mesh *mesh,
                                 const struct engine *e, struct gpart *gparts,

src/runner_doiact_grav.h

@@ -142,10 +142,9 @@ static INLINE void runner_dopair_grav_mm(struct runner *r,
 
   /* Some constants */
   const struct engine *e = r->e;
-  const struct space *s = e->s;
-  const int periodic = s->periodic;
-  const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
   const struct gravity_props *props = e->gravity_properties;
+  const int periodic = e->mesh->periodic;
+  const double dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]};
   const float r_s_inv = e->mesh->r_s_inv;
 
   TIMER_TIC;
@@ -681,12 +680,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
 
   /* Recover some useful constants */
   const struct engine *e = r->e;
-  const struct space *s = e->s;
-  const int periodic = s->periodic;
-  const float dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
-  const double r_s = e->mesh->r_s;
+  const int periodic = e->mesh->periodic;
+  const float dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]};
   const float r_s_inv = e->mesh->r_s_inv;
-  const double min_trunc = e->gravity_properties->r_cut_min * r_s;
+  const double min_trunc = e->mesh->r_cut_min;
 
   TIMER_TIC;
 
@@ -880,6 +877,7 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
 static INLINE void runner_doself_grav_pp_full(
     struct gravity_cache *restrict ci_cache, const int gcount,
     const int gcount_padded, const struct engine *e, struct gpart *gparts) {
+
   /* Loop over all particles in ci... */
   for (int pid = 0; pid < gcount; pid++) {
 
@@ -1089,11 +1087,9 @@ static INLINE void runner_doself_grav_pp(struct runner *r, struct cell *c) {
 
   /* Recover some useful constants */
   const struct engine *e = r->e;
-  const struct space *s = e->s;
-  const int periodic = s->periodic;
+  const int periodic = e->mesh->periodic;
   const float r_s_inv = e->mesh->r_s_inv;
-  const double r_s = e->mesh->r_s;
-  const double min_trunc = e->gravity_properties->r_cut_min * r_s;
+  const double min_trunc = e->mesh->r_cut_min;
 
   TIMER_TIC;
 
@@ -1185,12 +1181,11 @@ static INLINE void runner_dopair_recursive_grav(struct runner *r,
 
   /* Some constants */
   const struct engine *e = r->e;
-  const struct space *s = e->s;
   const int nodeID = e->nodeID;
-  const int periodic = s->periodic;
-  const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
+  const int periodic = e->mesh->periodic;
+  const double dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]};
   const double theta_crit2 = e->gravity_properties->theta_crit2;
-  const double max_distance = e->mesh->r_s * e->gravity_properties->r_cut_max;
+  const double max_distance = e->mesh->r_cut_max;
 
   /* Anything to do here? */
   if (!((cell_is_active_gravity(ci, e) && ci->nodeID == nodeID) ||
@@ -1391,24 +1386,18 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
 
   /* Some constants */
   const struct engine *e = r->e;
-  const struct space *s = e->s;
-  const int periodic = s->periodic;
-  const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
-  const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
+  const int periodic = e->mesh->periodic;
+  const double dim[3] = {e->mesh->dim[0], e->mesh->dim[1], e->mesh->dim[2]};
   const double theta_crit2 = e->gravity_properties->theta_crit2;
-  const double max_distance = e->mesh->r_s * e->gravity_properties->r_cut_max;
+  const double max_distance = e->mesh->r_cut_max;
 
-  /* const float u = 1.; */
-  /* float W; */
-  /* kernel_long_grav_force_eval(u, &W); */
-  /* message("a_smooth = %e, max_dist = %e, W=%e", */
-  /* 	  e->mesh->a_smooth, max_distance,W); */
+  const int cdim[3] = {e->s->cdim[0], e->s->cdim[1], e->s->cdim[2]};
 
   TIMER_TIC;
 
   /* Recover the list of top-level cells */
-  struct cell *cells = s->cells_top;
-  const int nr_cells = s->nr_cells;
+  struct cell *cells = e->s->cells_top;
+  const int nr_cells = e->s->nr_cells;
 
   /* Anything to do here? */
   if (!cell_is_active_gravity(ci, e)) return;