diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml
index f55ecc856953d4cb60a86e3461625318a1757693..346d2c0627ce2fdc1147d0d34fd4faab25b76559 100644
--- a/examples/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_6/eagle_6.yml
@@ -30,7 +30,7 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
- epsilon: 0.0001 # Softening length (in internal units).
+ 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 e619ddf246d46d4258dd66e6f9301eaf2dca7b94..f000b824ee134ad02c8dcf41b1ee12acf5ff7914 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1679,7 +1679,7 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
const struct engine *e = sp->e;
const int periodic = sp->periodic;
const double dim[3] = {sp->dim[0], sp->dim[1], sp->dim[2]};
- const double theta_crit = e->gravity_properties->theta_crit;
+ const double theta_crit2 = e->gravity_properties->theta_crit2;
/* Self interaction? */
if (cj == NULL) {
@@ -1733,7 +1733,7 @@ void cell_activate_subcell_grav_tasks(struct cell *ci, struct cell *cj,
const double r2 = dx * dx + dy * dy + dz * dz;
/* Can we use multipoles ? */
- if (gravity_multipole_accept(multi_i, multi_j, theta_crit, r2)) {
+ if (gravity_multipole_accept(multi_i, multi_j, theta_crit2, r2)) {
/* Ok, no need to drift anything */
return;
diff --git a/src/engine.c b/src/engine.c
index 8a66b56007811e2a1327f9b64fddfb478b9e650f..6f1864918875cec97214ffea60cb3bcfd5c64c46 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1726,7 +1726,7 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
const int cdim_ghost[3] = {s->cdim[0] / 4 + 1, s->cdim[1] / 4 + 1,
s->cdim[2] / 4 + 1};
- const double theta_crit = e->gravity_properties->theta_crit;
+ const double theta_crit2 = e->gravity_properties->theta_crit2;
struct cell *cells = s->cells_top;
const int n_ghosts = cdim_ghost[0] * cdim_ghost[1] * cdim_ghost[2] * 2;
@@ -1798,7 +1798,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 (!gravity_multipole_accept_rebuild(multi_i, multi_j, theta_crit,
+ if (!gravity_multipole_accept_rebuild(multi_i, multi_j, theta_crit2,
r2)) {
/* Ok, we need to add a direct pair calculation */
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index 18cf044434f7840a5a76f483540bb924a2365e26..8110b218eb064cdfc488afee9e0cdfb0c4f056f4 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -52,6 +52,8 @@ 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.");
+ p->theta_crit2 = p->theta_crit * p->theta_crit;
p->theta_crit_inv = 1. / p->theta_crit;
/* Softening lengths */
diff --git a/src/gravity_properties.h b/src/gravity_properties.h
index 2a5e4cb1e07ea591e2e3821704ec55abe7980360..c663f628b48c1af8d701792c256d7d30db248137 100644
--- a/src/gravity_properties.h
+++ b/src/gravity_properties.h
@@ -51,6 +51,9 @@ struct gravity_props {
/*! Tree opening angle (Multipole acceptance criterion) */
double theta_crit;
+ /*! Square of opening angle */
+ double theta_crit2;
+
/*! Inverse of opening angle */
double theta_crit_inv;
diff --git a/src/multipole.h b/src/multipole.h
index cf5b0ad9c90ad9c3abb773e53159d7c82f784afc..5d9211241483d7a1fa7789600eac6c5bcba9a45f 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -234,6 +234,7 @@ INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
* @brief Zeroes all the fields of a field tensor
*
* @param l The field tensor.
+ * @param ti_current The current (integer) time (for debugging only).
*/
INLINE static void gravity_field_tensors_init(struct grav_tensor *l,
integertime_t ti_current) {
@@ -246,7 +247,7 @@ INLINE static void gravity_field_tensors_init(struct grav_tensor *l,
}
/**
- * @brief Adds field tensrs to other ones (i.e. does la += lb).
+ * @brief Adds a field tensor to another one (i.e. does la += lb).
*
* @param la The gravity tensors to add to.
* @param lb The gravity tensors to add.
@@ -2174,6 +2175,9 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
l_b->F_102 += m_a->M_000 * D_soft_102(dx, dy, dz, r, eps_inv);
l_b->F_012 += m_a->M_000 * D_soft_012(dx, dy, dz, r, eps_inv);
l_b->F_111 += m_a->M_000 * D_soft_111(dx, dy, dz, r, eps_inv);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
+#error "Missing implementation for order >5"
#endif
}
}
@@ -2654,22 +2658,24 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
*
* @param ma The #multipole of the first #cell.
* @param mb The #multipole of the second #cell.
- * @param theta_crit The critical opening angle.
+ * @param theta_crit2 The square of the critical opening angle.
* @param r2 Square of the distance (periodically wrapped) between the
* multipoles.
*/
__attribute__((always_inline)) INLINE static int
gravity_multipole_accept_rebuild(const struct gravity_tensors *const ma,
const struct gravity_tensors *const mb,
- double theta_crit, double r2) {
+ double theta_crit2, double r2) {
- const double r_crit_a = ma->r_max_rebuild * theta_crit;
- const double r_crit_b = mb->r_max_rebuild * theta_crit;
+ const double r_crit_a = ma->r_max_rebuild;
+ const double r_crit_b = mb->r_max_rebuild;
+ const double size = r_crit_a + r_crit_b;
+ const double size2 = size * size;
// MATTHIEU: Make this mass-dependent ?
/* Multipole acceptance criterion (Dehnen 2002, eq.10) */
- return (r2 > (r_crit_a + r_crit_b) * (r_crit_a + r_crit_b));
+ return (r2 * theta_crit2 > size2);
}
/**
@@ -2681,21 +2687,23 @@ gravity_multipole_accept_rebuild(const struct gravity_tensors *const ma,
*
* @param ma The #multipole of the first #cell.
* @param mb The #multipole of the second #cell.
- * @param theta_crit The critical opening angle.
+ * @param theta_crit2 The square of the critical opening angle.
* @param r2 Square of the distance (periodically wrapped) between the
* multipoles.
*/
__attribute__((always_inline)) INLINE static int gravity_multipole_accept(
const struct gravity_tensors *const ma,
- const struct gravity_tensors *const mb, double theta_crit, double r2) {
+ const struct gravity_tensors *const mb, double theta_crit2, double r2) {
- const double r_crit_a = ma->r_max * theta_crit;
- const double r_crit_b = mb->r_max * theta_crit;
+ const double r_crit_a = ma->r_max;
+ const double r_crit_b = mb->r_max;
+ const double size = r_crit_a + r_crit_b;
+ const double size2 = size * size;
// MATTHIEU: Make this mass-dependent ?
/* Multipole acceptance criterion (Dehnen 2002, eq.10) */
- return (r2 > (r_crit_a + r_crit_b) * (r_crit_a + r_crit_b));
+ return (r2 * theta_crit2 > size2);
}
#endif /* SWIFT_MULTIPOLE_H */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 4c5f6a5f8b12da55a147c39f00825752c4ebf257..e4bfcc6e050750ff6b39e38d88e307cd1e7f22c9 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -1412,7 +1412,7 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
const double cell_width = s->width[0];
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const struct gravity_props *props = e->gravity_properties;
- const double theta_crit = props->theta_crit;
+ const double theta_crit2 = props->theta_crit2;
const double max_distance = props->a_smooth * props->r_cut_max * cell_width;
const double max_distance2 = max_distance * max_distance;
@@ -1473,7 +1473,7 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
* option... */
/* Can we use M-M interactions ? */
- if (gravity_multipole_accept(multi_i, multi_j, theta_crit, r2)) {
+ if (gravity_multipole_accept(multi_i, multi_j, theta_crit2, r2)) {
/* MATTHIEU: make a symmetric M-M interaction function ! */
runner_dopair_grav_mm(r, ci, cj);
@@ -1638,7 +1638,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
const int periodic = s->periodic;
const double cell_width = s->width[0];
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
- const double theta_crit = props->theta_crit;
+ const double theta_crit2 = props->theta_crit2;
const double max_distance = props->a_smooth * props->r_cut_max * cell_width;
const double max_distance2 = max_distance * max_distance;
@@ -1697,7 +1697,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
}
/* Check the multipole acceptance criterion */
- if (gravity_multipole_accept(multi_i, multi_j, theta_crit, r2)) {
+ if (gravity_multipole_accept(multi_i, multi_j, theta_crit2, r2)) {
/* Go for a (non-symmetric) M-M calculation */
runner_dopair_grav_mm(r, ci, cj);
@@ -1720,7 +1720,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
const double r2_rebuild = dx * dx + dy * dy + dz * dz;
/* Is the criterion violated now but was OK at the last rebuild ? */
- if (gravity_multipole_accept_rebuild(multi_i, multi_j, theta_crit,
+ if (gravity_multipole_accept_rebuild(multi_i, multi_j, theta_crit2,
r2_rebuild)) {
/* Alright, we have to take charge of that pair in a different way. */