diff --git a/src/cell.c b/src/cell.c
index d4cea2b09ef9dc5a102e46dd216a17dc34b604c4..a60347c4fca9218ded2d4c5275d8a20ca4bf1f6b 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1165,6 +1165,16 @@ void cell_check_multipole(struct cell *c, void *data) {
gravity_multipole_print(&c->multipole->m_pole);
error("Aborting");
}
+
+ /* Check that the upper limit of r_max is good enough */
+ if (!(c->multipole->r_max >= ma.r_max)) {
+ error("Upper-limit r_max=%e too small. Should be >=%e.",
+ c->multipole->r_max, ma.r_max);
+ } else if (c->multipole->r_max * c->multipole->r_max >
+ 3. * c->width[0] * c->width[0]) {
+ error("r_max=%e larger than cell diagonal %e.", c->multipole->r_max,
+ sqrt(3. * c->width[0] * c->width[0]));
+ }
}
#else
error("Calling debugging code without debugging flag activated.");
diff --git a/src/multipole.h b/src/multipole.h
index 60b00a459b807944ac1c736a02daa2b769e02e7c..c1d9d605e08aa044750c7e515d268c7f47c0631e 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -176,6 +176,9 @@ struct gravity_tensors {
/*! Centre of mass of the matter dsitribution */
double CoM[3];
+ /*! Upper limit of the CoM<->gpart distance */
+ double r_max;
+
/*! Multipole mass */
struct multipole m_pole;
@@ -967,9 +970,9 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
double com[3] = {0.0, 0.0, 0.0};
double vel[3] = {0.f, 0.f, 0.f};
- /* Collect the particle data. */
+ /* Collect the particle data for CoM. */
for (int k = 0; k < gcount; k++) {
- const float m = gparts[k].mass;
+ const double m = gparts[k].mass;
mass += m;
com[0] += gparts[k].x[0] * m;
@@ -989,7 +992,8 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
vel[1] *= imass;
vel[2] *= imass;
-/* Prepare some local counters */
+ /* Prepare some local counters */
+ double r_max2 = 0.;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
double M_100 = 0., M_010 = 0., M_001 = 0.;
#endif
@@ -1026,11 +1030,15 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
/* Construce the higher order terms */
for (int k = 0; k < gcount; k++) {
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
- const double m = gparts[k].mass;
const double dx[3] = {gparts[k].x[0] - com[0], gparts[k].x[1] - com[1],
gparts[k].x[2] - com[2]};
+ /* Maximal distance CoM<->gpart */
+ r_max2 = max(r_max2, dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ const double m = gparts[k].mass;
+
/* 1st order terms */
M_100 += -m * X_100(dx);
M_010 += -m * X_010(dx);
@@ -1115,6 +1123,7 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
/* Store the data on the multipole. */
m->m_pole.M_000 = mass;
+ m->r_max = sqrt(r_max2);
m->CoM[0] = com[0];
m->CoM[1] = com[1];
m->CoM[2] = com[2];
diff --git a/src/space.c b/src/space.c
index d9ce5f5a582a47cb79057c9e60787e2a0f64714b..e8e95c6ac53489907d3dbb22e0121f0c1a250587 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2083,15 +2083,38 @@ void space_split_recursive(struct space *s, struct cell *c,
/* Now shift progeny multipoles and add them up */
struct multipole temp;
+ double r_max = 0.;
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL) {
const struct cell *cp = c->progeny[k];
const struct multipole *m = &cp->multipole->m_pole;
+
+ /* Contribution to multipole */
gravity_M2M(&temp, m, c->multipole->CoM, cp->multipole->CoM,
s->periodic);
gravity_multipole_add(&c->multipole->m_pole, &temp);
+
+ /* Upper limit of max CoM<->gpart distance */
+ const float dx = c->multipole->CoM[0] - cp->multipole->CoM[0];
+ const float dy = c->multipole->CoM[1] - cp->multipole->CoM[1];
+ const float dz = c->multipole->CoM[2] - cp->multipole->CoM[2];
+ const float r2 = dx * dx + dy * dy + dz * dz;
+ r_max = max(r_max, cp->multipole->r_max + sqrtf(r2));
}
}
+ /* Alternative upper limit of max CoM<->gpart distance */
+ const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
+ ? c->multipole->CoM[0] - c->loc[0]
+ : c->loc[0] + c->width[0] - c->multipole->CoM[0];
+ const double dy = c->multipole->CoM[1] > c->loc[1] + c->width[1] / 2.
+ ? c->multipole->CoM[1] - c->loc[1]
+ : c->loc[1] + c->width[1] - c->multipole->CoM[1];
+ const double dz = c->multipole->CoM[2] > c->loc[2] + c->width[2] / 2.
+ ? c->multipole->CoM[2] - c->loc[2]
+ : c->loc[2] + c->width[2] - c->multipole->CoM[2];
+
+ /* Take minimum of both limits */
+ c->multipole->r_max = min(r_max, sqrt(dx * dx + dy * dy + dz * dz));
}
}
@@ -2151,13 +2174,24 @@ void space_split_recursive(struct space *s, struct cell *c,
/* Construct the multipole and the centre of mass*/
if (s->gravity) {
- if (gcount > 0)
+ if (gcount > 0) {
gravity_P2M(c->multipole, c->gparts, c->gcount);
- else {
+ const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
+ ? c->multipole->CoM[0] - c->loc[0]
+ : c->loc[0] + c->width[0] - c->multipole->CoM[0];
+ const double dy = c->multipole->CoM[1] > c->loc[1] + c->width[1] / 2.
+ ? c->multipole->CoM[1] - c->loc[1]
+ : c->loc[1] + c->width[1] - c->multipole->CoM[1];
+ const double dz = c->multipole->CoM[2] > c->loc[2] + c->width[2] / 2.
+ ? c->multipole->CoM[2] - c->loc[2]
+ : c->loc[2] + c->width[2] - c->multipole->CoM[2];
+ c->multipole->r_max = sqrt(dx * dx + dy * dy + dz * dz);
+ } else {
gravity_multipole_init(&c->multipole->m_pole);
c->multipole->CoM[0] = c->loc[0] + c->width[0] / 2.;
c->multipole->CoM[1] = c->loc[1] + c->width[1] / 2.;
c->multipole->CoM[2] = c->loc[2] + c->width[2] / 2.;
+ c->multipole->r_max = 0.;
}
}
}