Use a repartition based on resampling the old space nodeIDs.

Add safety by using an initial partition technique that should work anyway.

src/partition.c

@@ -927,7 +927,10 @@ static int check_complete(struct space *s, int verbose, int nregions) {
   int failed = 0;
   for (int i = 0; i < nregions; i++) present[i] = 0;
   for (int i = 0; i < s->nr_cells; i++) {
-    present[s->cells[i].nodeID]++;
+    if (s->cells[i].nodeID <= nregions)
+      present[s->cells[i].nodeID]++;
+    else
+      message("Bad nodeID: %d", s->cells[i].nodeID);
   }
   for (int i = 0; i < nregions; i++) {
     if (!present[i]) {
@@ -938,3 +941,53 @@ static int check_complete(struct space *s, int verbose, int nregions) {
   free(present);
   return (!failed);
 }
+
+
+/**
+ * @brief Partition a space of cells based on another space of cells.
+ *
+ * The two spaces are expected to be at different cell sizes, so what we'd
+ * like to do is assign the second space to geometrically closest nodes
+ * of the first, with the effect of minimizing particle movement when
+ * rebuilding the second space from the first.
+ *
+ * Since two spaces cannot exist simultaneously the old space is actually
+ * required in a decomposed state. These are the old cells sizes and counts
+ * per dimension, along with a list of the old nodeIDs. The old nodeIDs are
+ * indexed by the cellid (see cell_getid()), so should be stored that way.
+ *
+ * On exit the new space cells will have their nodeIDs assigned.
+ *
+ * @param oldh the cell dimensions of old space.
+ * @param oldcdim number of cells per dimension in old space.
+ * @param oldnodeIDs the nodeIDs of cells in the old space, indexed by old cellid.
+ * @param s the space to be partitioned.
+ *
+ * @return 1 if the new space contains nodeIDs from all nodes, 0 otherwise.
+ */
+int partition_space_to_space(double *oldh, double *oldcdim, int *oldnodeIDs,
+                             struct space *s) {
+
+  /* Loop over all the new cells. */
+  int nr_nodes = 0;
+  for (int i = 0; i < s->cdim[0]; i++) {
+    for (int j = 0; j < s->cdim[1]; j++) {
+      for (int k = 0; k < s->cdim[2]; k++) {
+
+        /* Scale indices to old cell space. */
+        int ii = rint(i * s->ih[0] * oldh[0]);
+        int jj = rint(j * s->ih[1] * oldh[1]);
+        int kk = rint(k * s->ih[2] * oldh[2]);
+
+        int cid = cell_getid(s->cdim, i, j, k);
+        int oldcid = cell_getid(oldcdim, ii, jj, kk);
+        s->cells[cid].nodeID = oldnodeIDs[oldcid];
+
+        if (oldnodeIDs[oldcid] > nr_nodes) nr_nodes = oldnodeIDs[oldcid];
+      }
+    }
+  }
+
+  /* Check we have all nodeIDs present in the resample. */
+  return check_complete(s, 1, nr_nodes + 1);
+}

src/partition.h

@@ -56,4 +56,6 @@ void partition_repartition(enum repartition_type reparttype, int nodeID,
 void partition_initial_partition(struct partition *initial_partition,
                                  int nodeID, int nr_nodes, struct space *s);
 
+int partition_space_to_space(double *oldh, double *oldcdim, int *oldnodeID,
+                             struct space *s);
 #endif /* SWIFT_PARTITION_H */

src/space.c

 /*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+*               2016 Peter W. Draper (p.w.draper@durham.ac.uk)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
@@ -205,9 +206,33 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
    * global partition is recomputed and the particles redistributed.
    * Be prepared to do that. */
 #ifdef WITH_MPI
-  int partition = 0;
-  if (cdim[0] < s->cdim[0] || cdim[1] < s->cdim[1] || cdim[2] < s->cdim[2])
-      partition = 1;
+  double oldh[3];
+  double oldcdim[3];
+  int *oldnodeIDs = NULL;
+  if (cdim[0] < s->cdim[0] || cdim[1] < s->cdim[1] || cdim[2] < s->cdim[2]) {
+
+    /* Capture state of current space. */
+    oldcdim[0] = s->cdim[0];
+    oldcdim[1] = s->cdim[1];
+    oldcdim[2] = s->cdim[2];
+    oldh[0] = s->h[0];
+    oldh[1] = s->h[1];
+    oldh[2] = s->h[2];
+
+    if ((oldnodeIDs = (int *)malloc(sizeof(int) * s->nr_cells)) == NULL)
+      error("Failed to allocate temporary nodeIDs.");
+
+    int cid = 0;
+    for (int i = 0; i < s->cdim[0]; i++) {
+      for (int j = 0; j < s->cdim[1]; j++) {
+        for (int k = 0; k < s->cdim[2]; k++) {
+          cid = cell_getid(oldcdim, i, j, k);
+          oldnodeIDs[cid] = s->cells[cid].nodeID;
+        }
+      }
+    }
+  }
+
 #endif
 
   /* Do we need to re-build the upper-level cells? */
@@ -268,29 +293,36 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
     fflush(stdout);
 
 #ifdef WITH_MPI
-    /* XXX create an engine_resplit() function to use the same method as the
-     * initial partition. Fake for now. */
-    if (partition) {
-        if (s->e->nodeID == 0)
-            message("cell dimensions have decreased. Recalculating the "
-                    "global partition.");
-
-        /* Change the global partitioning. */
-        int grid[3];
-        factor(s->e->nr_nodes, &grid[0], &grid[1]);
-        factor(s->e->nr_nodes / grid[1], &grid[0], &grid[2]);
-        factor(grid[0] * grid[1], &grid[1], &grid[0]);
-
-        /* Run through the cells and set their nodeID. */
-        int ind[3];
-        for (int k = 0; k < s->nr_cells; k++) {
-            c = &s->cells[k];
-            for (int j = 0; j < 3; j++) ind[j] = c->loc[j] / s->dim[j] * grid[j];
-            c->nodeID = ind[0] + grid[0] * (ind[1] + grid[1] * ind[2]);
-        }
+    if (oldnodeIDs != NULL) {
+      /* We have changed the top-level cell dimension, so need to redistribute
+       * cells around the nodes. We repartition using the old space node
+       * positions as a grid to resample. */
+      if (s->e->nodeID == 0)
+        message("basic cell dimensions have increased - recalculating the "
+                "global partition.");
+
+      if (!partition_space_to_space(oldh, oldcdim, oldnodeIDs, s) ) {
+
+        /* Failed, try another technique that requires no settings. */
+        message("Failed to get a new partition, trying less optimal method");
+        struct partition initial_partition;
+#ifdef HAVE_METIS
+        initial_partition.type = INITPART_METIS_NOWEIGHT;
+#else
+        initial_partition.type = INITPART_VECTORIZE;
+#endif
+        partition_initial_partition(&initial_partition, s->e->nodeID,
+                                    s->e->nr_nodes, s);
+      }
+
+      /* Re-distribute the particles to their new nodes. */
+      engine_redistribute(s->e);
+
+      /* Make the proxies. */
+      engine_makeproxies(s->e);
 
-        /* Make the proxies. */
-        engine_makeproxies(s->e);
+      /* Finished with these. */
+      free(oldnodeIDs);
     }
 #endif
   } /* re-build upper-level cells? */