diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml
index cf00a17694e9848825d96e1e5267961908dbbbb9..6f52835abf64d29ec3c0188722ce54ac9feb77f2 100644
--- a/examples/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_6/eagle_6.yml
@@ -13,6 +13,10 @@ TimeIntegration:
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
+
+#Scheduler:
+# tasks_per_cell: 8
+
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
diff --git a/src/cell.c b/src/cell.c
index 4e0d51c6c9c50a3a5b4463bde8139cba476f8021..cbc35a15c1ded28a09d7e9ef2d12479c1c068823 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1692,18 +1692,6 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
const double dim[3] = {sp->dim[0], sp->dim[1], sp->dim[2]};
const double theta_crit2 = e->gravity_properties->theta_crit2;
- /* Store the current dx_max and h_max values. */
- ci->multipole->r_max_old = ci->multipole->r_max;
- ci->multipole->CoM_old[0] = ci->multipole->CoM[0];
- ci->multipole->CoM_old[1] = ci->multipole->CoM[1];
- ci->multipole->CoM_old[2] = ci->multipole->CoM[2];
- if (cj != NULL) {
- cj->multipole->r_max_old = cj->multipole->r_max;
- cj->multipole->CoM_old[0] = cj->multipole->CoM[0];
- cj->multipole->CoM_old[1] = cj->multipole->CoM[1];
- cj->multipole->CoM_old[2] = cj->multipole->CoM[2];
- }
-
/* Self interaction? */
if (cj == NULL) {
@@ -1736,6 +1724,9 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
/* Anything to do here? */
if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
+ if (ci->ti_old_multipole < e->ti_current) cell_drift_multipole(ci, e);
+ if (cj->ti_old_multipole < e->ti_current) cell_drift_multipole(cj, e);
+
/* Recover the multipole information */
struct gravity_tensors *const multi_i = ci->multipole;
struct gravity_tensors *const multi_j = cj->multipole;
diff --git a/src/const.h b/src/const.h
index c8060a2be51468a791e192a65a74f1a4d9bc8e30..bfa0d611070ea133f238de18cb85292b5417a95d 100644
--- a/src/const.h
+++ b/src/const.h
@@ -113,4 +113,6 @@
#define SOURCETERMS_NONE
//#define SOURCETERMS_SN_FEEDBACK
+//#define ICHECK -30896
+
#endif /* SWIFT_CONST_H */
diff --git a/src/engine.c b/src/engine.c
index 47766a44545812c4d86970c8a09c9cbd970218a9..7a397e1a00b331946602b73dae9c5a11351b429b 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2971,7 +2971,7 @@ int engine_estimate_nr_tasks(struct engine *e) {
#endif
}
if (e->policy & engine_policy_self_gravity) {
- n1 += 24;
+ n1 += 32;
n2 += 1;
#ifdef WITH_MPI
n2 += 2;
@@ -3070,13 +3070,6 @@ void engine_rebuild(struct engine *e, int clean_h_values) {
/* Re-build the tasks. */
engine_maketasks(e);
- /* Run through the tasks and mark as skip or not. */
- if (engine_marktasks(e))
- error("engine_marktasks failed after space_rebuild.");
-
-/* Print the status of the system */
-// if (e->verbose) engine_print_task_counts(e);
-
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all cells have been drifted to the current time.
* That can include cells that have not
@@ -3085,6 +3078,13 @@ void engine_rebuild(struct engine *e, int clean_h_values) {
e->policy & engine_policy_self_gravity);
#endif
+ /* Run through the tasks and mark as skip or not. */
+ if (engine_marktasks(e))
+ error("engine_marktasks failed after space_rebuild.");
+
+ /* Print the status of the system */
+ // if (e->verbose) engine_print_task_counts(e);
+
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
@@ -3692,8 +3692,8 @@ void engine_step(struct engine *e) {
if (e->policy & engine_policy_reconstruct_mpoles)
engine_reconstruct_multipoles(e);
- // else
- // engine_drift_top_multipoles(e);
+ else
+ engine_drift_top_multipoles(e);
}
/* Print the number of active tasks ? */
@@ -3901,7 +3901,7 @@ void engine_do_drift_top_multipoles_mapper(void *map_data, int num_elements,
if (c != NULL && c->nodeID == e->nodeID) {
/* Drift the multipole at this level only */
- cell_drift_multipole(c, e);
+ if (c->ti_old_multipole != e->ti_current) cell_drift_multipole(c, e);
}
}
}
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
index f0d145647ab3f973f3c0ffc2f995ee01d534bc72..b83f8c73b4678145f2ecd6d94b136b5aadffccbd 100644
--- a/src/gravity/Default/gravity_debug.h
+++ b/src/gravity/Default/gravity_debug.h
@@ -26,6 +26,10 @@ __attribute__((always_inline)) INLINE static void gravity_debug_particle(
"x=[%.5e,%.5e,%.5e], v_full=[%.5e,%.5e,%.5e], a=[%.5e,%.5e,%.5e]\n",
p->mass, p->epsilon, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0],
p->v_full[1], p->v_full[2], p->a_grav[0], p->a_grav[1], p->a_grav[2]);
+#ifdef SWIFT_DEBUG_CHECKS
+ printf("num_interacted=%lld ti_drift=%lld ti_kick=%lld\n", p->num_interacted,
+ p->ti_drift, p->ti_kick);
+#endif
}
#endif /* SWIFT_DEFAULT_GRAVITY_DEBUG_H */
diff --git a/src/multipole.h b/src/multipole.h
index cdb675f413290690e4da1643ec1cb8c91a30a1ea..d842081814b6ef7b6f490a1ba47a0152dd84d5a1 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -188,18 +188,12 @@ struct gravity_tensors {
/*! Centre of mass of the matter dsitribution */
double CoM[3];
- /*! Centre of mass of the matter dsitribution at the last time-step */
- double CoM_old[3];
-
/*! Centre of mass of the matter dsitribution at the last rebuild */
double CoM_rebuild[3];
/*! Upper limit of the CoM<->gpart distance */
double r_max;
- /*! Upper limit of the CoM<->gpart distance at the last time-step */
- double r_max_old;
-
/*! Upper limit of the CoM<->gpart distance at the last rebuild */
double r_max_rebuild;
};
@@ -238,7 +232,7 @@ INLINE static void gravity_drift(struct gravity_tensors *m, double dt,
m->CoM[2] += dz;
/* Conservative change in maximal radius containing all gpart */
- m->r_max = m->r_max_rebuild + 0. * x_diff;
+ m->r_max = m->r_max_rebuild + x_diff;
}
/**
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 2844661049f156df489a71c33d161b01d08c084b..6eb9baad2035f9e22b6c84f34d70fdcf7d6f9165 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -155,7 +155,7 @@ void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
#endif
/* Do we need to drift the multipole ? */
- if (cj->ti_old_multipole != e->ti_current) cell_drift_multipole(cj, e);
+ if (cj->ti_old_multipole != e->ti_current) error("Undrifted multipole");
/* Let's interact at this level */
gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
@@ -469,9 +469,9 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
/* Do we need to drift the multipoles ? */
if (cj_active && ci->ti_old_multipole != e->ti_current)
- cell_drift_multipole(ci, e);
+ error("Un-drifted multipole");
if (ci_active && cj->ti_old_multipole != e->ti_current)
- cell_drift_multipole(cj, e);
+ error("Un-drifted multipole");
/* Centre of the cell pair */
const double loc[3] = {ci->loc[0], // + 0. * ci->width[0],
@@ -970,9 +970,9 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
struct gravity_tensors *const multi_j = cj->multipole;
/* Get the distance between the CoMs */
- double dx = multi_i->CoM_old[0] - multi_j->CoM_old[0];
- double dy = multi_i->CoM_old[1] - multi_j->CoM_old[1];
- double dz = multi_i->CoM_old[2] - multi_j->CoM_old[2];
+ 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) {
@@ -999,8 +999,7 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
* option... */
/* Can we use M-M interactions ? */
- if (gravity_M2L_accept(multi_i->r_max_old, multi_j->r_max_old, theta_crit2,
- r2)) {
+ if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2)) {
/* MATTHIEU: make a symmetric M-M interaction function ! */
runner_dopair_grav_mm(r, ci, cj);
@@ -1013,8 +1012,8 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
/* Alright, we'll have to split and recurse. */
else {
- const double ri_max = multi_i->r_max_old;
- const double rj_max = multi_j->r_max_old;
+ const double ri_max = multi_i->r_max;
+ const double rj_max = multi_j->r_max;
/* Split the larger of the two cells and start over again */
if (ri_max > rj_max) {
@@ -1170,8 +1169,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
/* Recover the local multipole */
struct gravity_tensors *const multi_i = ci->multipole;
- const double CoM_i[3] = {multi_i->CoM_old[0], multi_i->CoM_old[1],
- multi_i->CoM_old[2]};
+ const double CoM_i[3] = {multi_i->CoM[0], multi_i->CoM[1], multi_i->CoM[2]};
const double CoM_rebuild_i[3] = {multi_i->CoM_rebuild[0],
multi_i->CoM_rebuild[1],
multi_i->CoM_rebuild[2]};
@@ -1187,65 +1185,58 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
/* Avoid stupid cases */
if (ci == cj || cj->gcount == 0) continue;
- /* Get the distance between the CoMs */
- double dx = CoM_i[0] - multi_j->CoM_old[0];
- double dy = CoM_i[1] - multi_j->CoM_old[1];
- double dz = CoM_i[2] - multi_j->CoM_old[2];
+ /* Get the distance between the CoMs at the last rebuild*/
+ double dx_r = CoM_rebuild_i[0] - multi_j->CoM_rebuild[0];
+ double dy_r = CoM_rebuild_i[1] - multi_j->CoM_rebuild[1];
+ double dz_r = CoM_rebuild_i[2] - multi_j->CoM_rebuild[2];
/* Apply BC */
if (periodic) {
- dx = nearest(dx, dim[0]);
- dy = nearest(dy, dim[1]);
- dz = nearest(dz, dim[2]);
+ dx_r = nearest(dx_r, dim[0]);
+ dy_r = nearest(dy_r, dim[1]);
+ dz_r = nearest(dz_r, dim[2]);
}
- const double r2 = dx * dx + dy * dy + dz * dz;
+ const double r2_rebuild = dx_r * dx_r + dy_r * dy_r + dz_r * dz_r;
- /* Are we beyond the distance where the truncated forces are 0 ?*/
- if (periodic && r2 > max_distance2) {
+ /* Are we in charge of this cell pair? */
+ if (gravity_M2L_accept(multi_i->r_max_rebuild, multi_j->r_max_rebuild,
+ theta_crit2, r2_rebuild)) {
-#ifdef SWIFT_DEBUG_CHECKS
- /* Need to account for the interactions we missed */
- multi_i->pot.num_interacted += multi_j->m_pole.num_gpart;
-#endif
- continue;
- }
-
- /* Check the multipole acceptance criterion */
- if (gravity_M2L_accept(multi_i->r_max_old, multi_j->r_max_old, theta_crit2,
- r2)) {
-
- /* Go for a (non-symmetric) M-M calculation */
- runner_dopair_grav_mm(r, ci, cj);
-
- } else {
-
- /* Let's check whether we need to still operate on this pair */
-
- /* Get the distance between the CoMs at the last rebuild*/
- double dx_rebuild = CoM_rebuild_i[0] - multi_j->CoM_rebuild[0];
- double dy_rebuild = CoM_rebuild_i[1] - multi_j->CoM_rebuild[1];
- double dz_rebuild = CoM_rebuild_i[2] - multi_j->CoM_rebuild[2];
+ /* Let's compute the current distance between the cell pair*/
+ double dx = CoM_i[0] - multi_j->CoM[0];
+ double dy = CoM_i[1] - multi_j->CoM[1];
+ double dz = CoM_i[2] - multi_j->CoM[2];
/* Apply BC */
if (periodic) {
- dx_rebuild = nearest(dx_rebuild, dim[0]);
- dy_rebuild = nearest(dy_rebuild, dim[1]);
- dz_rebuild = nearest(dz_rebuild, dim[2]);
+ dx = nearest(dx, dim[0]);
+ dy = nearest(dy, dim[1]);
+ dz = nearest(dz, dim[2]);
}
- const double r2_rebuild = dx_rebuild * dx_rebuild +
- dy_rebuild * dy_rebuild +
- dz_rebuild * dz_rebuild;
+ const double r2 = dx * dx + dy * dy + dz * dz;
+
+ /* Are we beyond the distance where the truncated forces are 0 ?*/
+ if (periodic && r2 > max_distance2) {
- /* Is the criterion violated now but was OK at the last rebuild ? */
- if (gravity_M2L_accept(multi_i->r_max_rebuild, multi_j->r_max_rebuild,
- theta_crit2, r2_rebuild)) {
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Need to account for the interactions we missed */
+ multi_i->pot.num_interacted += multi_j->m_pole.num_gpart;
+#endif
+ continue;
+ }
+ /* Check the multipole acceptance criterion */
+ if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2)) {
+
+ /* Go for a (non-symmetric) M-M calculation */
+ runner_dopair_grav_mm(r, ci, cj);
+ } else {
/* Alright, we have to take charge of that pair in a different way. */
// MATTHIEU: We should actually open the tree-node here and recurse.
runner_dopair_grav_mm(r, ci, cj);
}
- }
- } /* Loop over top-level cells */
+ } /* We are in charge of this pair */
+ } /* Loop over top-level cells */
if (timer) TIMER_TOC(timer_dograv_long_range);
}