From 19e53963d71d527ee09156dd9e77d19e7d17ea0a Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sat, 10 Nov 2018 15:03:55 +0000
Subject: [PATCH] Self-consistent creation of gravity proxies and gravity
tasks.
---
src/engine.c | 45 ++++++++++++++++++++++++----------------
src/engine.h | 2 +-
src/engine_maketasks.c | 36 ++++++--------------------------
src/runner_doiact_grav.h | 34 +++++++++++++++---------------
4 files changed, 51 insertions(+), 66 deletions(-)
diff --git a/src/engine.c b/src/engine.c
index c125b09cd9..7f6b0ffd40 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3615,15 +3615,20 @@ void engine_makeproxies(struct engine *e) {
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] +
+ message("max_dist = %e", max_mesh_dist2);
+
+ /* Distance between centre of the cell and corners */
+ const double r_diag2 = cell_width[0] * cell_width[0] +
cell_width[1] * cell_width[1] +
cell_width[2] * cell_width[2];
-
- /* 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);
-
+ const double r_diag = 0.5 * sqrt(r_diag2);
+
+ /* Maximal distance from a shifted CoM to centre of cell */
+ const double delta_CoM = engine_max_proxy_centre_frac * r_diag;
+
+ /* Maximal distance from shifted CoM to any corner */
+ const double r_max = r_diag + 2. * delta_CoM;
+
/* Prepare the proxies and the proxy index. */
if (e->proxy_ind == NULL)
if ((e->proxy_ind = (int *)malloc(sizeof(int) * e->nr_nodes)) == NULL)
@@ -3632,20 +3637,20 @@ void engine_makeproxies(struct engine *e) {
e->nr_proxies = 0;
/* Compute how many cells away we need to walk */
- int delta = 1; /*hydro case */
+ int delta_cells = 1; /*hydro case */
/* Gravity needs to take the opening angle into account */
if (with_gravity) {
const double distance = 2. * r_max * theta_crit_inv;
- delta = (int)(distance / cells[0].dmin) + 1;
+ delta_cells = (int)(distance / cells[0].dmin) + 1;
}
/* Turn this into upper and lower bounds for loops */
- int delta_m = delta;
- int delta_p = delta;
+ int delta_m = delta_cells;
+ int delta_p = delta_cells;
/* Special case where every cell is in range of every other one */
- if (delta >= cdim[0] / 2) {
+ if (delta_cells >= cdim[0] / 2) {
if (cdim[0] % 2 == 0) {
delta_m = cdim[0] / 2;
delta_p = cdim[0] / 2 - 1;
@@ -3660,13 +3665,13 @@ void engine_makeproxies(struct engine *e) {
cells[i].num_grav_proxies = 0;
}
-
+ message("delta_m=%d delta_p=%d", delta_m, delta_p);
/* Let's be verbose about this choice */
if (e->verbose)
message(
"Looking for proxies up to %d top-level cells away (delta_m=%d "
"delta_p=%d)",
- delta, delta_m, delta_p);
+ delta_cells, delta_m, delta_p);
/* Loop over each cell in the space. */
for (int i = 0; i < cdim[0]; i++) {
@@ -3737,19 +3742,23 @@ void engine_makeproxies(struct engine *e) {
of cells cannot rely on just an M2L calculation. */
/* Minimal distance between any two points in the cells */
- const double min_dist2 = cell_min_dist2_same_size(&cells[cid], &cells[cjd], periodic, dim);
+ const double min_dist_centres2 = cell_min_dist2_same_size(&cells[cid], &cells[cjd], periodic, dim);
+
+ /* Let's now assume the CoMs will shift a bit */
+ const double min_dist_CoM = sqrt(min_dist_centres2) - 2. * delta_CoM;
+ const double min_dist_CoM2 = min_dist_CoM * min_dist_CoM;
/* 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_dist2 < max_mesh_dist2) &&
- (!gravity_M2L_accept(r_max, r_max, theta_crit2, min_dist2)))
+ if ((min_dist_CoM2 < max_mesh_dist2) &&
+ (!gravity_M2L_accept(r_max, r_max, theta_crit2, min_dist_CoM2)))
proxy_type |= (int)proxy_cell_type_gravity;
} else {
- if (!gravity_M2L_accept(r_max, r_max, theta_crit2, min_dist2))
+ if (!gravity_M2L_accept(r_max, r_max, theta_crit2, min_dist_CoM2))
proxy_type |= (int)proxy_cell_type_gravity;
}
}
diff --git a/src/engine.h b/src/engine.h
index 81e1b533b4..96d90a5e79 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -98,7 +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_max_proxy_centre_frac 0.2
#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 5daa1d1678..33ed4fa654 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -767,6 +767,7 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
struct cell *cells = s->cells_top;
const double theta_crit = e->gravity_properties->theta_crit;
const double max_distance = e->mesh->r_cut_max;
+ const double max_distance2 = max_distance * max_distance;
/* Compute how many cells away we need to walk */
const double distance = 2.5 * cells[0].width[0] / theta_crit;
@@ -807,17 +808,6 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci,
NULL);
-#ifdef SWIFT_DEBUG_CHECKS
- if (cell_getid(cdim, i, j, k) != cid)
- error("Incorrect calculation of indices (i,j,k)=(%d,%d,%d) cid=%d", i, j,
- k, cid);
-
- if (multi_i->r_max != multi_i->r_max_rebuild)
- error(
- "Multipole size not equal ot it's size after rebuild. But we just "
- "rebuilt...");
-#endif
-
/* Loop over every other cell within (Manhattan) range delta */
for (int ii = -delta_m; ii <= delta_p; ii++) {
int iii = i + ii;
@@ -836,7 +826,7 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
const int cjd = cell_getid(cdim, iii, jjj, kkk);
struct cell *cj = &cells[cjd];
- /* Avoid duplicates of local pairs*/
+ /* Avoid duplicates, empty cells and completely foreign pairs */
if (cid >= cjd || cj->grav.count == 0 ||
(ci->nodeID != nodeID && cj->nodeID != nodeID))
continue;
@@ -850,25 +840,11 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
if(multi_j == NULL && cj->nodeID != nodeID)
error("Multipole of cj was not exchanged properly via the proxies");
- /* Get the distance between the CoMs */
- double dx = multi_i->CoM[0] - multi_j->CoM[0];
- double dy = multi_i->CoM[1] - multi_j->CoM[1];
- double dz = multi_i->CoM[2] - multi_j->CoM[2];
-
- /* Apply BC */
- if (periodic) {
- dx = nearest(dx, dim[0]);
- dy = nearest(dy, dim[1]);
- dz = nearest(dz, dim[2]);
- }
- const double r2 = dx * dx + dy * dy + dz * dz;
-
/* Minimal distance between any pair of particles */
- const double min_radius =
- sqrt(r2) - (multi_i->r_max + multi_j->r_max);
+ const double min_radius2 = cell_min_dist2_same_size(ci, cj, periodic, dim);
/* Are we beyond the distance where the truncated forces are 0 ?*/
- if (periodic && min_radius > max_distance) continue;
+ if (periodic && min_radius2 > max_distance2) continue;
/* Are the cells too close for a MM interaction ? */
if (!cell_can_use_pair_mm_rebuild(ci, cj, e, s)) {
@@ -2066,8 +2042,8 @@ void engine_maketasks(struct engine *e) {
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);
+ 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
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 5bcd57d643..afe77eb92c 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -1714,7 +1714,7 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
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_cut_max;
+ const double max_distance2 = e->mesh->r_cut_max * e->mesh->r_cut_max;
TIMER_TIC;
@@ -1759,24 +1759,11 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
/* Skip empty cells */
if (multi_j->m_pole.M_000 == 0.f) continue;
- /* Get the distance between the CoMs at the last rebuild*/
- double dx_r = CoM_rebuild_top[0] - multi_j->CoM_rebuild[0];
- double dy_r = CoM_rebuild_top[1] - multi_j->CoM_rebuild[1];
- double dz_r = CoM_rebuild_top[2] - multi_j->CoM_rebuild[2];
-
- /* Apply BC */
- if (periodic) {
- dx_r = nearest(dx_r, dim[0]);
- dy_r = nearest(dy_r, dim[1]);
- dz_r = nearest(dz_r, dim[2]);
- }
- const double r2_rebuild = dx_r * dx_r + dy_r * dy_r + dz_r * dz_r;
-
- const double max_radius =
- sqrt(r2_rebuild) - (multi_top->r_max_rebuild + multi_j->r_max_rebuild);
+ /* Minimal distance between any pair of particles */
+ const double min_radius2 = cell_min_dist2_same_size(ci, cj, periodic, dim);
/* Are we beyond the distance where the truncated forces are 0 ?*/
- if (periodic && max_radius > max_distance) {
+ if (periodic && min_radius2 > max_distance2) {
#ifdef SWIFT_DEBUG_CHECKS
/* Need to account for the interactions we missed */
@@ -1790,6 +1777,19 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
continue;
}
+ /* Get the distance between the CoMs at the last rebuild*/
+ double dx_r = CoM_rebuild_top[0] - multi_j->CoM_rebuild[0];
+ double dy_r = CoM_rebuild_top[1] - multi_j->CoM_rebuild[1];
+ double dz_r = CoM_rebuild_top[2] - multi_j->CoM_rebuild[2];
+
+ /* Apply BC */
+ if (periodic) {
+ dx_r = nearest(dx_r, dim[0]);
+ dy_r = nearest(dy_r, dim[1]);
+ dz_r = nearest(dz_r, dim[2]);
+ }
+ const double r2_rebuild = dx_r * dx_r + dy_r * dy_r + dz_r * dz_r;
+
/* Are we in charge of this cell pair? */
if (gravity_M2L_accept(multi_top->r_max_rebuild, multi_j->r_max_rebuild,
theta_crit2, r2_rebuild)) {
--
GitLab