Drift the multipoles correctly and update r_max using the cell's x_diff

src/cell.c

@@ -2330,7 +2330,8 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
   if (ti_current < ti_old_multipole) error("Attempt to drift to the past");
 
   /* Drift the multipole */
-  if (ti_current > ti_old_multipole) gravity_drift(c->multipole, dt);
+  if (ti_current > ti_old_multipole)
+    gravity_drift(c->multipole, dt, c->dx_max_gpart);
 
   /* Are we not in a leaf ? */
   if (c->split) {
@@ -2365,7 +2366,8 @@ void cell_drift_multipole(struct cell *c, const struct engine *e) {
   /* Check that we are actually going to move forward. */
   if (ti_current < ti_old_multipole) error("Attempt to drift to the past");
 
-  if (ti_current > ti_old_multipole) gravity_drift(c->multipole, dt);
+  if (ti_current > ti_old_multipole)
+    gravity_drift(c->multipole, dt, c->dx_max_gpart);
 
   /* Update the time of the last drift */
   c->ti_old_multipole = ti_current;

src/multipole.h

@@ -216,8 +216,11 @@ INLINE static void gravity_reset(struct gravity_tensors *m) {
  *
  * @param m The #multipole.
  * @param dt The drift time-step.
+ * @param x_diff The maximal distance moved by any particle since the last
+ * rebuild.
  */
-INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
+INLINE static void gravity_drift(struct gravity_tensors *m, double dt,
+                                 float x_diff) {
 
   const double dx = m->m_pole.vel[0] * dt;
   const double dy = m->m_pole.vel[1] * dt;
@@ -229,8 +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 */
-  /* MATTHIEU: Use gpart->x_diff here ? */
-  m->r_max += sqrt(dx * dx + dy * dy + dz * dz);
+  m->r_max = m->r_max_rebuild + 2. * x_diff;
 }
 
 /**