Merge branch 'gravity_fixes' of gitlab.cosma.dur.ac.uk:swift/swiftsim into gravity_fixes

src/engine.c

@@ -4470,16 +4470,16 @@ void engine_makeproxies(struct engine *e) {
         /* Get the cell ID. */
         const int cid = cell_getid(cdim, i, j, k);
 
-	if (with_gravity) {
-
-	  /* Get ci's multipole */
-	  const struct gravity_tensors *multi_i = cells[cid].multipole;
-	  CoM_i[0] = multi_i->CoM[0];
-	  CoM_i[1] = multi_i->CoM[1];
-	  CoM_i[2] = multi_i->CoM[2];
-	  r_max_i = multi_i->r_max;
-	}
-	  
+        if (with_gravity) {
+
+          /* Get ci's multipole */
+          const struct gravity_tensors *multi_i = cells[cid].multipole;
+          CoM_i[0] = multi_i->CoM[0];
+          CoM_i[1] = multi_i->CoM[1];
+          CoM_i[2] = multi_i->CoM[2];
+          r_max_i = multi_i->r_max;
+        }
+
         /* Loop over every other cell */
         for (int ii = 0; ii < cdim[0]; ii++) {
           for (int jj = 0; jj < cdim[1]; jj++) {
@@ -4488,44 +4488,43 @@ void engine_makeproxies(struct engine *e) {
               /* Get the cell ID. */
               const int cjd = cell_getid(cdim, ii, jj, kk);
 
-	      /* In the gravity case, check distances using the MAC. */
-	      if (with_gravity) {
-		
-		/* Get cj's multipole */
-		const struct gravity_tensors *multi_j = cells[cjd].multipole;
-		const double CoM_j[3] = {multi_j->CoM[0], multi_j->CoM[1],
-					 multi_j->CoM[2]};
-		const double r_max_j = multi_j->r_max;
-		
-		/* Let's compute the current distance between the cell pair*/
-		double dx = CoM_i[0] - CoM_j[0];
-		double dy = CoM_i[1] - CoM_j[1];
-		double dz = CoM_i[2] - CoM_j[2];
-		
-		/* Apply BC */
-		if (periodic) {
-		  dx = nearest(dx, dim[0]);
-		  dy = nearest(dy, dim[1]);
-		  dz = nearest(dz, dim[2]);
-		}
-		const double r2 = dx * dx + dy * dy + dz * dz;
-
-		if (gravity_M2L_accept(r_max_i, r_max_j, theta_crit2, r2))
-		  continue;
-	      }
-
-	      /* In the hydro case, only carre about neighbours */
-	      else if(with_hydro) {
-
-		if (!((abs(i - ii) <= 1 || abs(i - ii - cdim[0]) <= 1 ||
-		     abs(i - ii + cdim[0]) <= 1) &&
-		    (abs(j - jj) <= 1 || abs(j - jj - cdim[1]) <= 1 ||
-		     abs(j - jj + cdim[1]) <= 1) &&
-		    (abs(k - kk) <= 1 || abs(k - kk - cdim[2]) <= 1 ||
-		     abs(k - kk + cdim[2]) <= 1))) 
-		  continue;
-	      }
+              /* In the gravity case, check distances using the MAC. */
+              if (with_gravity) {
+
+                /* Get cj's multipole */
+                const struct gravity_tensors *multi_j = cells[cjd].multipole;
+                const double CoM_j[3] = {multi_j->CoM[0], multi_j->CoM[1],
+                                         multi_j->CoM[2]};
+                const double r_max_j = multi_j->r_max;
+
+                /* Let's compute the current distance between the cell pair*/
+                double dx = CoM_i[0] - CoM_j[0];
+                double dy = CoM_i[1] - CoM_j[1];
+                double dz = CoM_i[2] - CoM_j[2];
+
+                /* Apply BC */
+                if (periodic) {
+                  dx = nearest(dx, dim[0]);
+                  dy = nearest(dy, dim[1]);
+                  dz = nearest(dz, dim[2]);
+                }
+                const double r2 = dx * dx + dy * dy + dz * dz;
 
+                if (gravity_M2L_accept(r_max_i, r_max_j, theta_crit2, r2))
+                  continue;
+              }
+
+              /* In the hydro case, only carre about neighbours */
+              else if (with_hydro) {
+
+                if (!((abs(i - ii) <= 1 || abs(i - ii - cdim[0]) <= 1 ||
+                       abs(i - ii + cdim[0]) <= 1) &&
+                      (abs(j - jj) <= 1 || abs(j - jj - cdim[1]) <= 1 ||
+                       abs(j - jj + cdim[1]) <= 1) &&
+                      (abs(k - kk) <= 1 || abs(k - kk - cdim[2]) <= 1 ||
+                       abs(k - kk + cdim[2]) <= 1)))
+                  continue;
+              }
 
               /* Add to proxies? */
               if (cells[cid].nodeID == e->nodeID &&

src/partition.c

@@ -354,23 +354,56 @@ static void pick_metis(struct space *s, int nregions, double *vertexw,
   /* Init the vertex weights array. */
   if (vertexw != NULL) {
     for (int k = 0; k < ncells; k++) {
-      if (vertexw[k] > 0) {
+      if (vertexw[k] > 1) {
         weights_v[k] = vertexw[k];
       } else {
         weights_v[k] = 1;
       }
     }
+
+#ifdef SWIFT_DEBUG_CHECKS
+    /* Check weights are all in range. */
+    int failed = 0;
+    for (int k = 0; k < ncells; k++) {
+      if ((idx_t)vertexw[k] < 0) {
+        message("Input vertex weight out of range: %ld", (long)vertexw[k]);
+        failed++;
+      }
+      if (weights_v[k] < 1) {
+        message("Used vertex weight  out of range: %" PRIDX, weights_v[k]);
+        failed++;
+      }
+    }
+    if (failed > 0) error("%d vertex weights are out of range", failed);
+#endif
   }
 
   /* Init the edges weights array. */
   if (edgew != NULL) {
     for (int k = 0; k < 26 * ncells; k++) {
-      if (edgew[k] > 0) {
+      if (edgew[k] > 1) {
         weights_e[k] = edgew[k];
       } else {
         weights_e[k] = 1;
       }
     }
+
+#ifdef SWIFT_DEBUG_CHECKS
+    /* Check weights are all in range. */
+    int failed = 0;
+    for (int k = 0; k < 26 * ncells; k++) {
+
+      if ((idx_t)edgew[k] < 0) {
+        message("Input edge weight out of range: %ld", (long)edgew[k]);
+        failed++;
+      }
+      if (weights_e[k] < 1) {
+        message("Used edge weight out of range: %" PRIDX, weights_e[k]);
+        failed++;
+      }
+    }
+    if (failed > 0) error("%d edge weights are out of range", failed);
+#endif
   }
 
   /* Set the METIS options. */
@@ -661,35 +694,42 @@ static void repart_edge_metis(int partweights, int bothweights, int timebins,
 
     double wmine = weights_e[0];
     double wmaxe = weights_e[0];
-    for (int k = 0; 26 * k < nr_cells; k++) {
+    for (int k = 0; k < 26 * nr_cells; k++) {
       wmaxe = weights_e[k] > wmaxe ? weights_e[k] : wmaxe;
       wmine = weights_e[k] < wmine ? weights_e[k] : wmine;
     }
 
-    double wscalev = 1.0;
-    double wscalee = 1.0;
     if (bothweights) {
 
-      /* Make maximum value same in both weights systems. */
-      if (wmaxv > wmaxe) {
-        wscalee = wmaxv / wmaxe;
+      /* Make range the same in both weights systems. */
+      if ((wmaxv - wminv) > (wmaxe - wmine)) {
+        double wscale = (wmaxv - wminv) / (wmaxe - wmine);
+        for (int k = 0; k < 26 * nr_cells; k++) {
+          weights_e[k] = (weights_e[k] - wmine) * wscale + wminv;
+        }
+        wmine = wminv;
         wmaxe = wmaxv;
+
       } else {
-        wscalev = wmaxe / wmaxv;
+        double wscale = (wmaxe - wmine) / (wmaxv - wminv);
+        for (int k = 0; k < nr_cells; k++) {
+          weights_v[k] = (weights_v[k] - wminv) * wscale + wmine;
+        }
+        wminv = wmine;
         wmaxv = wmaxe;
       }
 
       /* Scale to the METIS range. */
-      wscalev *= metis_maxweight / (wmaxv - wminv);
+      double wscale = (metis_maxweight - 1.0) / (wmaxv - wminv);
       for (int k = 0; k < nr_cells; k++) {
-        weights_v[k] = (weights_v[k] - wminv) * wscalev + 1.0;
+        weights_v[k] = (weights_v[k] - wminv) * wscale + 1.0;
       }
     }
 
     /* Scale to the METIS range. */
-    wscalee *= metis_maxweight / (wmaxe - wmine);
+    double wscale = (metis_maxweight - 1.0) / (wmaxe - wmine);
     for (int k = 0; k < 26 * nr_cells; k++) {
-      weights_e[k] = (weights_e[k] - wmine) * wscalee + 1.0;
+      weights_e[k] = (weights_e[k] - wmine) * wscale + 1.0;
     }
 
     /* And partition, use both weights or not as requested. */