diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml
index 3bfe4b53f5582625cf0ba1a4eefed13cac2f9c71..cf00a17694e9848825d96e1e5267961908dbbbb9 100644
--- a/examples/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_6/eagle_6.yml
@@ -13,9 +13,6 @@ 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:
- cell_split_size: 64
-
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
@@ -29,8 +26,8 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
- theta: 0.7 # Opening angle (Multipole acceptance criterion)
- epsilon: 0.001 # Softening length (in internal units).
+ theta: 0.85 # Opening angle (Multipole acceptance criterion)
+ epsilon: 0.001 # Softening length (in internal units).
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/src/cell.c b/src/cell.c
index 692caa89e29aac2f27c6bafcaf6d8d9687e2d6aa..8e9fd3b11d25c057e62ef854a8020e49f76dd0e0 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1692,6 +1692,19 @@ 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) {
diff --git a/src/engine.c b/src/engine.c
index 47766a44545812c4d86970c8a09c9cbd970218a9..6de5ae303e66eaa2afaad0e46a18b81e459b0541 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3656,7 +3656,7 @@ void engine_step(struct engine *e) {
if (e->nodeID == 0) {
/* Print some information to the screen */
- printf(" %6d %14e %14e %10zu %10zu %10zu %21.3f\n", e->step, e->time,
+ printf(" %6d %10lld %3d %14e %14e %10zu %10zu %10zu %21.3f\n", e->step, e->ti_current, e->max_active_bin, e->time,
e->timeStep, e->updates, e->g_updates, e->s_updates,
e->wallclock_time);
fflush(stdout);
diff --git a/src/gravity_cache.h b/src/gravity_cache.h
index fdc89605765b460b355b3958e34287991be5ff1b..097e4e1cb36e74ddb548fb3ab74196dd66512e6e 100644
--- a/src/gravity_cache.h
+++ b/src/gravity_cache.h
@@ -189,13 +189,14 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
* that have a reasonable magnitude. We use the cell width for this */
const float pos_padded[3] = {-2. * cell->width[0], -2. * cell->width[1],
-2. * cell->width[2]};
+ const float eps_padded = epsilon[0];
/* Pad the caches */
for (int i = gcount; i < gcount_padded; ++i) {
x[i] = pos_padded[0];
y[i] = pos_padded[1];
z[i] = pos_padded[2];
- epsilon[i] = 0.f;
+ epsilon[i] = eps_padded;
m[i] = 0.f;
active[i] = 0;
use_mpole[i] = 0;
@@ -248,12 +249,14 @@ gravity_cache_populate_no_mpole(timebin_t max_active_bin,
* that have a reasonable magnitude. We use the cell width for this */
const float pos_padded[3] = {-2. * cell->width[0], -2. * cell->width[1],
-2. * cell->width[2]};
+ const float eps_padded = epsilon[0];
+
/* Pad the caches */
for (int i = gcount; i < gcount_padded; ++i) {
x[i] = pos_padded[0];
y[i] = pos_padded[1];
z[i] = pos_padded[2];
- epsilon[i] = 0.f;
+ epsilon[i] = eps_padded;
m[i] = 0.f;
active[i] = 0;
}
diff --git a/src/multipole.h b/src/multipole.h
index 370b7d52d11f040fa053e45dd69ee4f2416eb0b2..cdb675f413290690e4da1643ec1cb8c91a30a1ea 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -188,12 +188,18 @@ 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;
};
@@ -232,7 +238,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 + 2. * x_diff;
+ m->r_max = m->r_max_rebuild + 0. * x_diff;
}
/**
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index ca814ed4bc29a1a79571627704635e2994ffc220..ed68476a625fb8eed14fee38bba438a11863b301 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -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[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];
+ 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];
/* Apply BC */
if (periodic) {
@@ -999,7 +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, multi_j->r_max, theta_crit2, r2)) {
+ if (gravity_M2L_accept(multi_i->r_max_old, multi_j->r_max_old, theta_crit2, r2)) {
/* MATTHIEU: make a symmetric M-M interaction function ! */
runner_dopair_grav_mm(r, ci, cj);
@@ -1012,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;
- const double rj_max = multi_j->r_max;
+ const double ri_max = multi_i->r_max_old;
+ const double rj_max = multi_j->r_max_old;
/* Split the larger of the two cells and start over again */
if (ri_max > rj_max) {
@@ -1169,7 +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[0], multi_i->CoM[1], multi_i->CoM[2]};
+ const double CoM_i[3] = {multi_i->CoM_old[0], multi_i->CoM_old[1], multi_i->CoM_old[2]};
const double CoM_rebuild_i[3] = {multi_i->CoM_rebuild[0],
multi_i->CoM_rebuild[1],
multi_i->CoM_rebuild[2]};
@@ -1186,9 +1186,9 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
if (ci == cj || cj->gcount == 0) continue;
/* Get the distance between the CoMs */
- 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];
+ 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];
/* Apply BC */
if (periodic) {
@@ -1209,7 +1209,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
}
/* Check the multipole acceptance criterion */
- if (gravity_M2L_accept(multi_i->r_max, multi_j->r_max, theta_crit2, r2)) {
+ 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);
diff --git a/src/vector.h b/src/vector.h
index e463bbd2e62f0c6b51446d22a9805670b134e684..3222a8a9b6a1bf7f5dc02e914a168d7c5c911b68 100644
--- a/src/vector.h
+++ b/src/vector.h
@@ -433,7 +433,7 @@ __attribute__((always_inline)) INLINE vector vec_reciprocal_sqrt(vector x) {
#else
/* Needed for cache alignment. */
-#define VEC_SIZE 16
+#define VEC_SIZE 8
#endif /* WITH_VECTORIZATION */
#endif /* VEC_MACRO */