Most methods now working, MPI implementations added

examples/main.c

@@ -81,9 +81,17 @@ int main(int argc, char *argv[]) {
   char ICfileName[200] = "";
   float dt_max = 0.0f;
   ticks tic;
-  int nr_nodes = 1, myrank = 0, grid[3] = {1, 1, 1};
+  int nr_nodes = 1, myrank = 0;
   FILE *file_thread;
   int with_outputs = 1;
+  struct pgrid pgrid;
+
+  /* Default parition type is grid. */
+  pgrid.type = GRID_GRID;
+  pgrid.grid[0] = 1;
+  pgrid.grid[1] = 1;
+  pgrid.grid[2] = 1;
+
 
 /* Choke on FP-exceptions. */
 // feenableexcept( FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW );
@@ -107,9 +115,9 @@ int main(int argc, char *argv[]) {
   fflush(stdout);
 
   /* Set a default grid so that grid[0]*grid[1]*grid[2] == nr_nodes. */
-  factor(nr_nodes, &grid[0], &grid[1]);
-  factor(nr_nodes / grid[1], &grid[0], &grid[2]);
-  factor(grid[0] * grid[1], &grid[1], &grid[0]);
+  factor(nr_nodes, &pgrid.grid[0], &pgrid.grid[1]);
+  factor(nr_nodes / pgrid.grid[1], &pgrid.grid[0], &pgrid.grid[2]);
+  factor(pgrid.grid[0] * pgrid.grid[1], &pgrid.grid[1], &pgrid.grid[0]);
 #endif
 
   /* Greeting message */
@@ -142,8 +150,30 @@ int main(int argc, char *argv[]) {
         if (!strcpy(ICfileName, optarg)) error("Error parsing IC file name.");
         break;
       case 'g':
-        if (sscanf(optarg, "%i %i %i", &grid[0], &grid[1], &grid[2]) != 3)
-          error("Error parsing grid.");
+        /* Grid is one of "g", "r", "m", "w", or "v". g can be followed by three
+         * numbers. */
+        switch (optarg[0]) {
+          case 'g':
+            pgrid.type = GRID_GRID;
+            if (strlen(optarg) > 2) {
+              if (sscanf(optarg, "g %i %i %i", &pgrid.grid[0], &pgrid.grid[1],
+                         &pgrid.grid[2]) != 3) 
+                error("Error parsing grid.");
+            }
+            break;
+          case 'r':
+            pgrid.type = GRID_RANDOM;
+            break;
+          case 'm':
+            pgrid.type = GRID_METIS_NOWEIGHT;
+            break;
+          case 'w':
+            pgrid.type = GRID_METIS_WEIGHT;
+            break;
+          case 'v':
+            pgrid.type = GRID_VECTORIZE;
+            break;
+        }
         break;
       case 'm':
         if (sscanf(optarg, "%lf", &h_max) != 1) error("Error parsing h_max.");
@@ -193,7 +223,7 @@ int main(int argc, char *argv[]) {
 #if defined(WITH_MPI)
   if (myrank == 0) {
     message("Running with %i thread(s) per node.", nr_threads);
-    message("grid set to [ %i %i %i ].", grid[0], grid[1], grid[2]);
+    message("grid set to [ %i %i %i ].", pgrid.grid[0], pgrid.grid[1], pgrid.grid[2]);
 
     if (nr_nodes == 1) {
       message("WARNING: you are running with one MPI rank.");
@@ -326,7 +356,7 @@ int main(int argc, char *argv[]) {
 
 #ifdef WITH_MPI
   /* Split the space. */
-  engine_split(&e, grid);
+  engine_split(&e, &pgrid);
   engine_redistribute(&e);
 #endif
 
@@ -380,7 +410,8 @@ int main(int argc, char *argv[]) {
 
 /* Repartition the space amongst the nodes? */
 #if defined(WITH_MPI) && defined(HAVE_METIS)
-    if (j % 100 == 2) e.forcerepart = 1;
+    //if (j % 100 == 2) e.forcerepart = 1;
+    if (j % 10 == 9) e.forcerepart = 1;
 #endif
 
     timers_reset(timers_mask_all);

src/Makefile.am

@@ -35,13 +35,13 @@ endif
 # List required headers
 include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
-    common_io.h single_io.h multipole.h map.h tools.h poisson_disc.h
+    common_io.h single_io.h multipole.h map.h tools.h partition.h
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
     serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
     units.c common_io.c single_io.c multipole.c version.c map.c \
-    kernel.c tools.c poisson_disc.c
+    kernel.c tools.c partition.c
 
 # Include files for distribution, not installation.
 noinst_HEADERS = atomic.h cycle.h error.h inline.h kernel.h vector.h \

src/engine.c

@@ -49,7 +49,7 @@
 #include "debug.h"
 #include "error.h"
 #include "timers.h"
-#include "poisson_disc.h"
+#include "partition.h"
 
 #ifdef LEGACY_GADGET2_SPH
 #include "runner_iact_legacy.h"
@@ -2067,52 +2067,142 @@ void engine_makeproxies(struct engine *e) {
 /**
  * @brief Split the underlying space according to the given grid.
  *
+ * Only used with MPI the partitions produced associate cells with nodes.
+ *
  * @param e The #engine.
  * @param grid The grid.
  */
-
-void engine_split(struct engine *e, int *grid) {
+void engine_split(struct engine *e, struct pgrid *grid) {
 
 #ifdef WITH_MPI
 
-  //int j, k;
-  //int ind[3];
   struct space *s = e->s;
-  //struct cell *c;
-
-  /* If we've got the wrong number of nodes, fail. */
-  //if (e->nr_nodes != grid[0] * grid[1] * grid[2])
-  //   error("Grid size does not match number of nodes.");
-
-  /* Run through the cells and set their nodeID. */
-  // message("s->dim = [%e,%e,%e]", s->dim[0], s->dim[1], s->dim[2]);
-  //for (k = 0; k < s->nr_cells; k++) {
-  //  c = &s->cells[k];
-  //  for (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]);
-  //  // message("cell at [%e,%e,%e]: ind = [%i,%i,%i], nodeID = %i", c->loc[0],
-  //  // c->loc[1], c->loc[2], ind[0], ind[1], ind[2], c->nodeID);
-  //}
-
-  /* Poisson split is stochastic so can only be done by one node. */
-  float *samplelist = malloc(sizeof( float ) * e->nr_nodes * 3);
-  if ( samplelist == NULL )
-    error("Failed to allocate samplelist");
-
-  if (e->nodeID == 0) {
-    if ( poisson_generate(s, e->nr_nodes, samplelist) == 0 )
-      error("Failed to partition cells");
-    message("samplelist[0,1,2] = %f,%f,%f", samplelist[0], samplelist[1], samplelist[2]);
+
+  if (grid->type == GRID_GRID) {
+    int j, k;
+    int ind[3];
+    struct cell *c;
+
+    /* If we've got the wrong number of nodes, fail. */
+    if (e->nr_nodes != grid->grid[0] * grid->grid[1] * grid->grid[2])
+      error("Grid size does not match number of nodes.");
+
+    /* Run through the cells and set their nodeID. */
+    // message("s->dim = [%e,%e,%e]", s->dim[0], s->dim[1], s->dim[2]);
+    for (k = 0; k < s->nr_cells; k++) {
+      c = &s->cells[k];
+      for (j = 0; j < 3; j++) ind[j] = c->loc[j] / s->dim[j] * grid->grid[j];
+      c->nodeID = ind[0] + grid->grid[0] * (ind[1] + grid->grid[1] * ind[2]);
+      // message("cell at [%e,%e,%e]: ind = [%i,%i,%i], nodeID = %i", c->loc[0],
+      // c->loc[1], c->loc[2], ind[0], ind[1], ind[2], c->nodeID);
+    }
   }
+  else if (grid->type == GRID_RANDOM) {
+
+    /* Random can only be done by one node, each node requires a list of
+     * the sampling positions. */
+    float *samplelist = malloc(sizeof( float ) * e->nr_nodes * 4);
+    if ( samplelist == NULL )
+      error("Failed to allocate samplelist");
 
-  /* Share the samplelist of points around all the nodes. */
-  int res = MPI_Bcast(samplelist, e->nr_nodes * 3, MPI_FLOAT, 0, MPI_COMM_WORLD);
-  if (res != MPI_SUCCESS)
-    mpi_error(res,"Failed to bcast the partition samples.");
+    if (e->nodeID == 0) {
+      if ( part_pick_random(s, e->nr_nodes, samplelist) == 0 )
+        error("Failed to partition cells");
+    }
+
+    /* Share the samplelist of points around all the nodes. */
+    int res = MPI_Bcast(samplelist, e->nr_nodes * 4, MPI_FLOAT, 0,
+                        MPI_COMM_WORLD);
+    if (res != MPI_SUCCESS)
+      mpi_error(res,"Failed to bcast the partition samples.");
+
+    /* And apply to our cells */
+    part_split_random(s, e->nr_nodes, samplelist);
+    free(samplelist);
+  }
+  else if (grid->type == GRID_METIS_WEIGHT || grid->type == GRID_METIS_NOWEIGHT) {
+
+    /* Simple k-way partition selected by METIS using cell particle counts as
+     * weights or not. Should be best when starting with a inhomogeneous dist. */
+
+    /* Space for particles per cell counts, which will be used as weights or not. */
+    int *weights = NULL;
+    if (grid->type == GRID_METIS_WEIGHT) {
+      if ((weights = malloc(sizeof(int) * s->nr_cells)) == NULL)
+        error("Failed to allocate weights buffer.");
+      bzero(weights, sizeof(int) * s->nr_cells);
+
+      /* Check each particle and accumilate the counts per cell. */
+      struct part *parts = s->parts;
+      int *cdim = s->cdim;
+      double ih[3], dim[3];
+      ih[0] = s->ih[0];
+      ih[1] = s->ih[1];
+      ih[2] = s->ih[2];
+      dim[0] = s->dim[0];
+      dim[1] = s->dim[1];
+      dim[2] = s->dim[2];
+      for (int k = 0; k < s->nr_parts; k++) {
+        for (int j = 0; j < 3; j++) {
+          if (parts[k].x[j] < 0.0)
+            parts[k].x[j] += dim[j];
+          else if (parts[k].x[j] >= dim[j])
+            parts[k].x[j] -= dim[j];
+        }
+        const int cid = cell_getid(cdim, parts[k].x[0] * ih[0],
+                                   parts[k].x[1] * ih[1], parts[k].x[2] * ih[2]);
+        weights[cid]++;
+      }
+      
+      /* Get all the counts from all the nodes. */
+      if (MPI_Allreduce(MPI_IN_PLACE, weights, s->nr_cells, MPI_INT, MPI_SUM,
+                        MPI_COMM_WORLD) != MPI_SUCCESS)
+        error("Failed to allreduce particle cell weights.");
+    }
 
-  /* And apply to our cells */
-  poisson_split(s, e->nr_nodes, samplelist);
-  free(samplelist);
+
+    /* Main node does the partition calculation. */
+    int *celllist = malloc(sizeof( int ) * s->nr_cells);
+    if (celllist == NULL)
+      error("Failed to allocate celllist");
+
+    if (e->nodeID == 0)
+      part_pick_metis(s, e->nr_nodes, weights, celllist);
+
+    /* Distribute the celllist partition and apply. */
+    int res = MPI_Bcast(celllist, s->nr_cells, MPI_INT, 0, MPI_COMM_WORLD);
+    if (res != MPI_SUCCESS)
+      mpi_error(res,"Failed to bcast the cell list");
+
+    /* And apply to our cells */
+    part_split_metis(s, e->nr_nodes, celllist);
+
+    if (weights != NULL) free(weights);
+    free(celllist);
+  }
+  else if (grid->type == GRID_VECTORIZE) {
+
+    /* Vectorised selection, guaranteed to work, but not very clumpy in the
+     * selection of regions. */
+    int *samplecells = malloc(sizeof( int ) * e->nr_nodes * 3);
+    if ( samplecells == NULL )
+      error("Failed to allocate samplecells");
+
+    if (e->nodeID == 0) {
+      part_pick_vector(s, e->nr_nodes, samplecells);
+    }
+
+    /* Share the samplecells around all the nodes. */
+    int res = MPI_Bcast(samplecells, e->nr_nodes * 3, MPI_INT, 0,
+                        MPI_COMM_WORLD);
+    if (res != MPI_SUCCESS)
+      mpi_error(res,"Failed to bcast the partition sample cells.");
+
+    /* And apply to our cells */
+    part_split_vector(s, e->nr_nodes, samplecells);
+    free(samplecells);
+
+  }
 
   /* Make the proxies. */
   engine_makeproxies(e);
@@ -2136,7 +2226,7 @@ void engine_split(struct engine *e, int *grid) {
   s->parts = parts_new;
   s->xparts = xparts_new;
 
-#endif /* WITH_MPI */
+#endif
 }
 
 /**

src/engine.h

@@ -29,6 +29,7 @@
 #include "scheduler.h"
 #include "space.h"
 #include "task.h"
+#include "param.h"
 
 /* Some constants. */
 #define engine_policy_none 0
@@ -135,7 +136,7 @@ void engine_launch(struct engine *e, int nr_runners, unsigned int mask);
 void engine_prepare(struct engine *e);
 void engine_step(struct engine *e);
 void engine_maketasks(struct engine *e);
-void engine_split(struct engine *e, int *grid);
+void engine_split(struct engine *e, struct pgrid *grid);
 int engine_exchange_strays(struct engine *e, int offset, int *ind, int N);
 void engine_rebuild(struct engine *e);
 void engine_repartition(struct engine *e);

src/param.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2015 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
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PARAM_H
+#define SWIFT_PARAM_H
+
+
+/* Initial partition grid struct. Defines type of partitioning to use and any
+ * related parameters. */
+enum grid_types {
+  GRID_GRID = 0,
+  GRID_RANDOM,
+  GRID_VECTORIZE,
+  GRID_METIS_WEIGHT,
+  GRID_METIS_NOWEIGHT
+};
+
+struct pgrid {
+  enum grid_types type;
+  int grid[3];
+};
+
+
+#endif /* SWIFT_PARAM_H */

src/partition.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2015 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
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PARTITION_H
+#define SWIFT_PARTITION_H
+
+#include "space.h"
+#include "cell.h"
+
+int part_pick_random(struct space *s, int nregions, float *samplelist);
+void part_split_random(struct space *s, int nregions, float *samplelist);
+
+void part_pick_vector(struct space *s, int nregions, int *samplecells);
+void part_split_vector(struct space *s, int nregions, int *samplecells);
+
+void part_pick_metis(struct space *s, int nregions, int *weight, int *celllist);
+void part_split_metis(struct space *s, int nregions, int *celllist);
+
+#endif /* SWIFT_POISSON_DISC_H */

src/swift.h

@@ -46,6 +46,8 @@
 #include "timers.h"
 #include "units.h"
 #include "tools.h"
+#include "param.h"
+#include "partition.h"
 #include "version.h"
 
 #ifdef LEGACY_GADGET2_SPH