diff --git a/src/cell.h b/src/cell.h
index 59780697f42bce4b1fb778d7fcdf471de9cb2b27..7fda5fb40c69e1bc83a367265a0fbe9afb9ff266 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -670,6 +670,73 @@ int cell_can_use_pair_mm(const struct cell *ci, const struct cell *cj,
int cell_can_use_pair_mm_rebuild(const struct cell *ci, const struct cell *cj,
const struct engine *e, const struct space *s);
+/**
+ * @brief Compute the sqaure of the minimal distance between any two points in two cells of the same size
+ *
+ * @param ci The first #cell.
+ * @param cj The second #cell.
+ * @param periodic Are we using periodic BCs?
+ * @param dim The dimensions of the simulation volume
+ */
+__attribute__((always_inline)) INLINE static double
+cell_min_dist2_same_size(const struct cell *restrict ci, const struct cell *restrict cj, const int periodic,
+ const double dim[3]) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci->width[0] != cj->width[0])
+ error("Cells of different size!");
+ if (ci->width[1] != cj->width[1])
+ error("Cells of different size!");
+ if (ci->width[2] != cj->width[2])
+ error("Cells of different size!");
+#endif
+
+ const double cix_min = ci->loc[0];
+ const double ciy_min = ci->loc[1];
+ const double ciz_min = ci->loc[2];
+ const double cjx_min = cj->loc[0];
+ const double cjy_min = cj->loc[1];
+ const double cjz_min = cj->loc[2];
+
+ const double cix_max = ci->loc[0] + ci->width[0];
+ const double ciy_max = ci->loc[1] + ci->width[1];
+ const double ciz_max = ci->loc[2] + ci->width[2];
+ const double cjx_max = cj->loc[0] + cj->width[0];
+ const double cjy_max = cj->loc[1] + cj->width[1];
+ const double cjz_max = cj->loc[2] + cj->width[2];
+
+ if (periodic) {
+
+ const double dx = min4(fabs(nearest(cix_min - cjx_min, dim[0])),
+ fabs(nearest(cix_min - cjx_max, dim[0])),
+ fabs(nearest(cix_max - cjx_min, dim[0])),
+ fabs(nearest(cix_max - cjx_max, dim[0])));
+
+ const double dy = min4(fabs(nearest(ciy_min - cjy_min, dim[1])),
+ fabs(nearest(ciy_min - cjy_max, dim[1])),
+ fabs(nearest(ciy_max - cjy_min, dim[1])),
+ fabs(nearest(ciy_max - cjy_max, dim[1])));
+
+ const double dz = min4(fabs(nearest(ciz_min - cjz_min, dim[2])),
+ fabs(nearest(ciz_min - cjz_max, dim[2])),
+ fabs(nearest(ciz_max - cjz_min, dim[2])),
+ fabs(nearest(ciz_max - cjz_max, dim[2])));
+
+ return dx * dx + dy * dy + dz * dz;
+
+ } else {
+
+ const double dx = min(fabs(cix_max - cjx_min),
+ fabs(cix_min - cjx_max));
+ const double dy = min(fabs(ciy_max - cjy_min),
+ fabs(ciy_min - cjy_max));
+ const double dz = min(fabs(ciz_max - cjz_min),
+ fabs(ciz_min - cjz_max));
+
+ return dx * dx + dy * dy + dz * dz;
+ }
+}
+
/* Inlined functions (for speed). */
/**
diff --git a/src/engine.c b/src/engine.c
index 5fddc4a319898bbfb3bc5a3dd53d85f1c4d9db7f..c125b09cd9b2a32a8192bfa833ad1fbd70f4579d 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3612,15 +3612,17 @@ void engine_makeproxies(struct engine *e) {
const int with_gravity = (e->policy & engine_policy_self_gravity);
const double theta_crit_inv = e->gravity_properties->theta_crit_inv;
const double theta_crit2 = e->gravity_properties->theta_crit2;
- const double max_distance = e->mesh->r_cut_max;
+ const double max_mesh_dist = e->mesh->r_cut_max;
+ const double max_mesh_dist2 = max_mesh_dist * max_mesh_dist;
/* Maximal distance between CoMs and any particle in the cell */
const double r_max2 = cell_width[0] * cell_width[0] +
cell_width[1] * cell_width[1] +
cell_width[2] * cell_width[2];
- const double r_max = sqrt(r_max2);
- message("PROXIES: r_max= %e max_distance= %e", r_max, max_distance);
+ /* Maximal distance assuming the CoM can have moved by from the centre by
+ engine_max_proxy_centre_frac (default: 20%) of the cell */
+ const double r_max = sqrt(r_max2) * 0.5 * (1. + 2. * engine_max_proxy_centre_frac);
/* Prepare the proxies and the proxy index. */
if (e->proxy_ind == NULL)
@@ -3673,11 +3675,7 @@ void engine_makeproxies(struct engine *e) {
/* Get the cell ID. */
const int cid = cell_getid(cdim, i, j, k);
-
- /* and it's location */
- const double loc_i[3] = {cells[cid].loc[0], cells[cid].loc[1],
- cells[cid].loc[2]};
-
+
/* Loop over all its neighbours neighbours in range. */
for (int ii = -delta_m; ii <= delta_p; ii++) {
int iii = i + ii;
@@ -3738,44 +3736,20 @@ void engine_makeproxies(struct engine *e) {
for an M2L interaction and hence require a proxy as this pair
of cells cannot rely on just an M2L calculation. */
- const double loc_j[3] = {cells[cjd].loc[0], cells[cjd].loc[1],
- cells[cjd].loc[2]};
-
- /* Start with the distance between the cell centres. */
- double dx = loc_i[0] - loc_j[0];
- double dy = loc_i[1] - loc_j[1];
- double dz = loc_i[2] - loc_j[2];
-
- /* Apply BC */
- if (periodic) {
- dx = nearest(dx, dim[0]);
- dy = nearest(dy, dim[1]);
- dz = nearest(dz, dim[2]);
- }
-
- /* Add to it for the case where the future CoMs are in the
- * corners */
- dx += cell_width[0];
- dy += cell_width[1];
- dz += cell_width[2];
-
- /* This is a crazy upper-bound but the best we can do */
- const double r2 = dx * dx + dy * dy + dz * dz;
-
- /* Minimal distance between any pair of particles */
- const double min_radius = sqrt(r2) - 2. * r_max;
+ /* Minimal distance between any two points in the cells */
+ const double min_dist2 = cell_min_dist2_same_size(&cells[cid], &cells[cjd], periodic, dim);
/* Are we beyond the distance where the truncated forces are 0
* but not too far such that M2L can be used? */
if (periodic) {
- if ((min_radius < max_distance) &&
- (!gravity_M2L_accept(r_max, r_max, theta_crit2, r2)))
+ if ((min_dist2 < max_mesh_dist2) &&
+ (!gravity_M2L_accept(r_max, r_max, theta_crit2, min_dist2)))
proxy_type |= (int)proxy_cell_type_gravity;
} else {
- if (!gravity_M2L_accept(r_max, r_max, theta_crit2, r2))
+ if (!gravity_M2L_accept(r_max, r_max, theta_crit2, min_dist2))
proxy_type |= (int)proxy_cell_type_gravity;
}
}
diff --git a/src/engine.h b/src/engine.h
index 50eef4b314da4c294a86e59e29542859f7f402f2..81e1b533b43f6c26214d30fc7c19976f00e3ea2f 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -98,6 +98,7 @@ enum engine_step_properties {
#define engine_maxproxies 64
#define engine_tasksreweight 1
#define engine_parts_size_grow 1.05
+#define engine_max_proxy_centre_frac 0.25
#define engine_redistribute_alloc_margin 1.2
#define engine_default_energy_file_name "energy"
#define engine_default_timesteps_file_name "timesteps"
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index c831b31ed6bcf3c237f945b1f65aa4100a6f2eaa..5daa1d167827dee5995f3e8218c1ac914be961e3 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -890,11 +890,13 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
const struct proxy *p = &e->proxies[proxy_id];
/* Check whether the cell exists in the proxy */
- int n = 0;
- for (n = 0; n < p->nr_cells_in; n++)
- if(p->cells_in[n] == cj)
+ int n = 0, err = 1;
+ for (; n < p->nr_cells_in; n++)
+ if(p->cells_in[n] == cj) {
+ err = 0;
break;
- if(n == p->nr_cells_in)
+ }
+ if(err)
error("Cell %d not found in the proxy but trying to construct grav task!",
cjd);
}
@@ -908,11 +910,13 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
const struct proxy *p = &e->proxies[proxy_id];
/* Check whether the cell exists in the proxy */
- int n = 0;
- for (n = 0; n < p->nr_cells_in; n++)
- if(p->cells_in[n] == ci)
+ int n = 0, err = 1;
+ for (; n < p->nr_cells_in; n++)
+ if(p->cells_in[n] == ci) {
+ err = 0;
break;
- if(n == p->nr_cells_in)
+ }
+ if(err)
error("Cell %d not found in the proxy but trying to construct grav task!",
cid);
}
@@ -2059,6 +2063,13 @@ void engine_maketasks(struct engine *e) {
message("Linking stars tasks took %.3f %s (including reweight).",
clocks_from_ticks(getticks() - tic2), clocks_getunit());
+ if(e->nodeID == 0)
+ for(int i = 0; i < e->s->nr_cells; ++i) {
+ message("cid= %d num_grav_proxy= %d num_hydro_proxy= %d num_foreign_grav= %d num_foreign_hydro= %d",
+ i, e->s->cells_top[i].num_grav_proxies, e->s->cells_top[i].num_hydro_proxies,
+ e->s->cells_top[i].num_foreign_pair_grav, e->s->cells_top[i].num_foreign_pair_hydro);
+ }
+
#ifdef WITH_MPI
if (e->policy & engine_policy_feedback)
error("Cannot run stellar feedback with MPI (yet).");