Also drift foreign top-level multipoles every time-step

src/engine.c

@@ -1977,8 +1977,8 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
           const int cjd = cell_getid(cdim, iii, jjj, kkk);
           struct cell *cj = &cells[cjd];
 
-	  if(i==11 && j==0 && k==10)
-	    message("Found direct neighbour: (i,j,k)=(%d,%d,%d) (iii,jjj,kkk)=(%d,%d,%d) nodeID=%d", 	i,j,k, iii,jjj,kkk, cj->nodeID);
+	  /* if(i==11 && j==0 && k==10) */
+	  /*   message("Found direct neighbour: (i,j,k)=(%d,%d,%d) (iii,jjj,kkk)=(%d,%d,%d) nodeID=%d", 	i,j,k, iii,jjj,kkk, cj->nodeID); */
 
 
 	  /* Avoid duplicates of local pairs*/
@@ -3947,11 +3947,11 @@ void engine_step(struct engine *e) {
   if (e->policy & engine_policy_self_gravity) {
     for (int i = 0; i < e->s->nr_cells; ++i)
       num_gpart_mpole += e->s->cells_top[i].multipole->m_pole.num_gpart;
-    if (num_gpart_mpole != e->s->nr_gparts)
+    if (num_gpart_mpole != e->total_nr_gparts)
       error(
           "Multipoles don't contain the total number of gpart mpoles=%zd "
           "ngparts=%zd",
-          num_gpart_mpole, e->s->nr_gparts);
+          num_gpart_mpole, e->total_nr_gparts);
   }
 #endif
 
@@ -4140,7 +4140,7 @@ void engine_do_drift_top_multipoles_mapper(void *map_data, int num_elements,
 
   for (int ind = 0; ind < num_elements; ind++) {
     struct cell *c = &cells[ind];
-    if (c != NULL && c->nodeID == e->nodeID) {
+    if (c != NULL) {
 
       /* Drift the multipole at this level only */
       if (c->ti_old_multipole != e->ti_current) cell_drift_multipole(c, e);

src/runner_doiact_grav.h

@@ -1241,8 +1241,8 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
       if(check) {
 	++direct_ngbs;
 	direct_ngbs_gpart += cj->multipole->m_pole.num_gpart;
-	message("Found direct neighbour %d: (i,j,k)=(%d,%d,%d) (ii,jj,kk)=(%d,%d,%d) nodeID=%d",
-		direct_ngbs, i,j,k, ii,jj,kk, cj->nodeID);
+	/* message("Found direct neighbour %d: (i,j,k)=(%d,%d,%d) (ii,jj,kk)=(%d,%d,%d) nodeID=%d", */
+	/* 	direct_ngbs, i,j,k, ii,jj,kk, cj->nodeID); */
       }
 #endif
       

src/space.c

@@ -234,6 +234,8 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->xparts = NULL;
     c->gparts = NULL;
     c->sparts = NULL;
+    if(s->gravity)
+      bzero(c->multipole, sizeof(struct gravity_tensors));
     for (int i = 0; i < 13; i++)
       if (c->sort[i] != NULL) {
         free(c->sort[i]);
@@ -264,6 +266,7 @@ void space_free_cells(struct space *s) {
   threadpool_map(&s->e->threadpool, space_rebuild_recycle_mapper, s->cells_top,
                  s->nr_cells, sizeof(struct cell), 0, s);
   s->maxdepth = 0;
+  message("Done"); fflush(stdout);
 }
 
 /**
@@ -953,6 +956,8 @@ void space_rebuild(struct space *s, int verbose) {
   // message( "hooking up cells took %.3f %s." ,
   // clocks_from_ticks(getticks() - tic), clocks_getunit());
 
+  message("ti_old=%lld", ti_old);
+
   /* At this point, we have the upper-level cells, old or new. Now make
      sure that the parts in each cell are ok. */
   space_split(s, cells_top, s->nr_cells, verbose);