diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml index ba181c8374d5aad935c257ca46b6997111fe2048..53eb99765d859a26dcdeb8c9148fb74c058c1544 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 0308442571bb76e8697bec10ceaa224022b96a91..f7e9a07c40c1620770355a9b6351a72ff515d259 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 9c7a31acc905b1c1321f489deb8fac27be74da88..dd7e1e267246c0f73e760191a071c6e24b96cfe8 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 1e2bf22e8f3d8db1fb49264704afca65449b6a81..892517d11c94ae42b8253e20dd958bf10bff92cb 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 43a1d0d78f403ba9d2f1202db5e8d7e648ad6a11..27a5de0a4102cae4ca787c10c60cf3bbc3a983ee 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 c6689874baf54de8e5c449282b5a8ddc04bee11d..e878d53290a4d62dd976c80dc77af8fd605e46c2 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 26d6b90c745e91554c87e5c476b228f50d14cb38..65089a23d57ed6c64352b3e22c2e201a97ad13c4 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 352056cf340e7bbbce336e03e67a060f172155b1..0f38ff3e9421c122b61caf74290c170b62b2dd92 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))