From cf17832d56159f9c4c80df642b00f9e4cc415744 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Thu, 23 Nov 2017 19:21:54 +0000
Subject: [PATCH] Removed more ddebugging code. Restored the gravity-pair
construction loop based on the MAC.
---
examples/EAGLE_12/eagle_12.yml | 2 +-
examples/EAGLE_6/eagle_6.yml | 2 +-
src/const.h | 3 --
src/engine.c | 68 ++++++++--------------------------
src/gravity_properties.c | 2 +-
src/runner.c | 6 +--
src/runner_doiact_grav.h | 2 +
src/timeline.h | 2 +-
8 files changed, 23 insertions(+), 64 deletions(-)
diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml
index ba181c8374..53eb99765d 100644
--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -27,7 +27,7 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
epsilon: 0.001 # Softening length (in internal units).
- theta: 2.5 # Opening angle (Multipole acceptance criterion)
+ theta: 0.85 # Opening angle (Multipole acceptance criterion)
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml
index 0308442571..f7e9a07c40 100644
--- a/examples/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_6/eagle_6.yml
@@ -30,7 +30,7 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
- theta: 2.5 # Opening angle (Multipole acceptance criterion)
+ theta: 0.85 # Opening angle (Multipole acceptance criterion)
epsilon: 0.001 # Softening length (in internal units).
# Parameters for the hydrodynamics scheme
diff --git a/src/const.h b/src/const.h
index 9c7a31acc9..dd7e1e2672 100644
--- a/src/const.h
+++ b/src/const.h
@@ -110,7 +110,4 @@
#define SOURCETERMS_NONE
//#define SOURCETERMS_SN_FEEDBACK
-//#define ICHECK 5726454604296ll
-#define ICHECK 67457606141961ll
-
#endif /* SWIFT_CONST_H */
diff --git a/src/engine.c b/src/engine.c
index 1e2bf22e8f..892517d11c 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2125,38 +2125,15 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
struct gravity_tensors *const multi_i = ci->multipole;
const double CoM_i[3] = {multi_i->CoM[0], multi_i->CoM[1], multi_i->CoM[2]};
- /* /\* Loop over every other cell *\/ */
- /* for (int ii = 0; ii < cdim[0]; ii++) { */
- /* for (int jj = 0; jj < cdim[1]; jj++) { */
- /* for (int kk = 0; kk < cdim[2]; kk++) { */
+ /* Loop over every other cell */
+ for (int ii = 0; ii < cdim[0]; ii++) {
+ for (int jj = 0; jj < cdim[1]; jj++) {
+ for (int kk = 0; kk < cdim[2]; kk++) {
- /* /\* Get the cell *\/ */
- /* const int cjd = cell_getid(cdim, ii, jj, kk); */
- /* struct cell *cj = &cells[cjd]; */
-
- /* Now loop over all the neighbours of this cell */
- for (int ii = -1; ii < 2; ii++) {
- int iii = i + ii;
- if (!s->periodic && (iii < 0 || iii >= cdim[0])) continue;
- iii = (iii + cdim[0]) % cdim[0];
- for (int jj = -1; jj < 2; jj++) {
- int jjj = j + jj;
- if (!s->periodic && (jjj < 0 || jjj >= cdim[1])) continue;
- jjj = (jjj + cdim[1]) % cdim[1];
- for (int kk = -1; kk < 2; kk++) {
- int kkk = k + kk;
- if (!s->periodic && (kkk < 0 || kkk >= cdim[2])) continue;
- kkk = (kkk + cdim[2]) % cdim[2];
-
- /* Get the neighbouring cell */
- const int cjd = cell_getid(cdim, iii, jjj, kkk);
+ /* Get the cell */
+ const int cjd = cell_getid(cdim, ii, jj, kk);
struct cell *cj = &cells[cjd];
- /* if(i==11 && j==0 && k==10) */
- /* message("Found direct neighbour: (i,j,k)=(%d,%d,%d)
- * (iii,jjj,kkk)=(%d,%d,%d) nodeID=%d", i,j,k, iii,jjj,kkk,
- * cj->nodeID); */
-
/* Avoid duplicates of local pairs*/
if (cid <= cjd && cj->nodeID == nodeID) continue;
@@ -2180,8 +2157,7 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
const double r2 = dx * dx + dy * dy + dz * dz;
/* Are the cells too close for a MM interaction ? */
- if (1 ||
- !gravity_M2L_accept(multi_i->r_max_rebuild,
+ if (!gravity_M2L_accept(multi_i->r_max_rebuild,
multi_j->r_max_rebuild, theta_crit2, r2)) {
/* Ok, we need to add a direct pair calculation */
@@ -3262,8 +3238,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
/* Periodic gravity stuff (Note this is not linked to a cell) ? */
- else if (t->type == task_type_grav_top_level) { //||
- // t->type == task_type_grav_ghost) {
+ else if (t->type == task_type_grav_top_level) {
scheduler_activate(s, t);
}
@@ -4484,6 +4459,8 @@ void engine_makeproxies(struct engine *e) {
struct cell *cells = s->cells_top;
struct proxy *proxies = e->proxies;
ticks tic = getticks();
+ //const int with_hydro = (e->policy & engine_policy_hydro);
+ //const int with_gravity = (e->policy & engine_policy_self_gravity);
/* Prepare the proxies and the proxy index. */
if (e->proxy_ind == NULL)
@@ -4505,25 +4482,12 @@ void engine_makeproxies(struct engine *e) {
/* Get the cell ID. */
const int cid = cell_getid(cdim, ind[0], ind[1], ind[2]);
- /* Loop over all its neighbours (periodic). */
- for (int i = -1; i <= 1; i++) {
- int ii = ind[0] + i;
- if (ii >= cdim[0])
- ii -= cdim[0];
- else if (ii < 0)
- ii += cdim[0];
- for (int j = -1; j <= 1; j++) {
- int jj = ind[1] + j;
- if (jj >= cdim[1])
- jj -= cdim[1];
- else if (jj < 0)
- jj += cdim[1];
- for (int k = -1; k <= 1; k++) {
- int kk = ind[2] + k;
- if (kk >= cdim[2])
- kk -= cdim[2];
- else if (kk < 0)
- kk += cdim[2];
+ /* Loop over every other cell */
+ for (int ii = 0; ii < cdim[0]; ii++) {
+ for (int jj = 0; jj < cdim[1]; jj++) {
+ for (int kk = 0; kk < cdim[2]; kk++) {
+
+ //if(with_hydro)
/* Get the cell ID. */
const int cjd = cell_getid(cdim, ii, jj, kk);
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index 43a1d0d78f..27a5de0a41 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -52,7 +52,7 @@ void gravity_props_init(struct gravity_props *p,
/* Opening angle */
p->theta_crit = parser_get_param_double(params, "Gravity:theta");
- // if (p->theta_crit >= 1.) error("Theta too large. FMM won't converge.");
+ if (p->theta_crit >= 1.) error("Theta too large. FMM won't converge.");
p->theta_crit2 = p->theta_crit * p->theta_crit;
p->theta_crit_inv = 1. / p->theta_crit;
diff --git a/src/runner.c b/src/runner.c
index c6689874ba..e878d53290 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -552,9 +552,6 @@ void runner_do_init_grav(struct runner *r, struct cell *c, int timer) {
/* Anything to do here? */
if (!cell_is_active(c, e)) return;
- /* Drift the multipole */
- // cell_drift_multipole(c, e);
-
/* Reset the gravity acceleration tensors */
gravity_field_tensors_init(&c->multipole->pot, e->ti_current);
@@ -1489,8 +1486,7 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
/* Check that this gpart has interacted with all the other
* particles (via direct or multipoles) in the box */
- if (gp->num_interacted != e->total_nr_gparts &&
- gp->id_or_neg_offset == ICHECK)
+ if (gp->num_interacted != e->total_nr_gparts)
error(
"g-particle (id=%lld, type=%s) did not interact "
"gravitationally "
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 26d6b90c74..65089a23d5 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -1301,12 +1301,14 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
error("Not found the right number of particles in top-level interactions");
#endif
+#ifdef ICHECK
if (check)
message(
"Interacted with %d indirectly and ignored %d direct interactions "
"(counter=%lld) nr_cells=%d total=%lld",
other_ngbs_gpart, direct_ngbs_gpart, counter, nr_cells,
e->total_nr_gparts);
+#endif
if (timer) TIMER_TOC(timer_dograv_long_range);
}
diff --git a/src/timeline.h b/src/timeline.h
index 352056cf34..0f38ff3e94 100644
--- a/src/timeline.h
+++ b/src/timeline.h
@@ -32,7 +32,7 @@ typedef long long integertime_t;
typedef char timebin_t;
/*! The number of time bins */
-#define num_time_bins 26
+#define num_time_bins 56
/*! The maximal number of timesteps in a simulation */
#define max_nr_timesteps (1LL << (num_time_bins + 1))
--
GitLab