Do not do the L2L or L2P if a field tensor has received 0 contributions.

src/multipole.h

@@ -100,6 +100,9 @@ struct grav_tensor {
   integertime_t ti_init;
 
 #endif
+
+  /* Has this tensor received any contribution? */
+  char interacted;
 };
 
 struct multipole {
@@ -216,18 +219,18 @@ INLINE static void gravity_reset(struct gravity_tensors *m) {
  */
 INLINE static void gravity_drift(struct gravity_tensors *m, double dt) {
 
-  /* const double dx = m->m_pole.vel[0] * dt; */
-  /* const double dy = m->m_pole.vel[1] * dt; */
-  /* const double dz = m->m_pole.vel[2] * dt; */
+  const double dx = m->m_pole.vel[0] * dt;
+  const double dy = m->m_pole.vel[1] * dt;
+  const double dz = m->m_pole.vel[2] * dt;
 
   /* Move the whole thing according to bulk motion */
-  /* m->CoM[0] += dx; */
-  /* m->CoM[1] += dy; */
-  /* m->CoM[2] += dz; */
+  m->CoM[0] += dx;
+  m->CoM[1] += dy;
+  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 += sqrt(dx * dx + dy * dy + dz * dz);
 }
 
 /**
@@ -259,6 +262,8 @@ INLINE static void gravity_field_tensors_add(struct grav_tensor *la,
   la->num_interacted += lb->num_interacted;
 #endif
 
+  la->interacted = 1;
+
   /* Add 0th order terms */
   la->F_000 += lb->F_000;
 
@@ -347,6 +352,7 @@ INLINE static void gravity_field_tensors_add(struct grav_tensor *la,
 INLINE static void gravity_field_tensors_print(const struct grav_tensor *l) {
 
   printf("-------------------------\n");
+  printf("Interacted: %d\n", l->interacted);
   printf("F_000= %12.5e\n", l->F_000);
 #if SELF_GRAVITY_MULTIPOLE_ORDER > 0
   printf("-------------------------\n");
@@ -1544,6 +1550,9 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
   l_b->num_interacted += m_a->num_gpart;
 #endif
 
+  /* Record that this tensor has received contributions */
+  l_b->interacted = 1;
+
   /*  0th order term */
   l_b->F_000 += m_a->M_000 * pot.D_000;
 

src/runner.c

@@ -1472,9 +1472,11 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
 #ifdef SWIFT_DEBUG_CHECKS
         if (e->policy & engine_policy_self_gravity) {
 
+          /* Let's add a self interaction to simplify the count */
+          gp->num_interacted++;
+
           /* Check that this gpart has interacted with all the other
            * particles (via direct or multipoles) in the box */
-          gp->num_interacted++;
           if (gp->num_interacted != (long long)e->s->nr_gparts)
             error(
                 "g-particle (id=%lld, type=%s) did not interact "

src/runner_doiact_grav.h

@@ -68,12 +68,16 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
 #endif
         struct grav_tensor shifted_tensor;
 
-        /* Shift the field tensor */
-        gravity_L2L(&shifted_tensor, &c->multipole->pot, cp->multipole->CoM,
-                    c->multipole->CoM);
+        /* If the tensor received any contribution, push it down */
+        if (c->multipole->pot.interacted) {
 
-        /* Add it to this level's tensor */
-        gravity_field_tensors_add(&cp->multipole->pot, &shifted_tensor);
+          /* Shift the field tensor */
+          gravity_L2L(&shifted_tensor, &c->multipole->pot, cp->multipole->CoM,
+                      c->multipole->CoM);
+
+          /* Add it to this level's tensor */
+          gravity_field_tensors_add(&cp->multipole->pot, &shifted_tensor);
+        }
 
         /* Recurse */
         runner_do_grav_down(r, cp, 0);
@@ -82,6 +86,9 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
 
   } else { /* Leaf case */
 
+    /* We can abort early if no interactions via multipole happened */
+    if (!c->multipole->pot.interacted) return;
+
     if (!cell_are_gpart_drifted(c, e)) error("Un-drifted gparts");
 
     /* Apply accelerations to the particles */
@@ -97,6 +104,8 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
         /* Check that particles have been drifted to the current time */
         if (gp->ti_drift != e->ti_current)
           error("gpart not drifted to current time");
+        if (c->multipole->pot.ti_init != e->ti_current)
+          error("c->field tensor not initialised");
 #endif
 
         /* Apply the kernel */