Merged with latest master.

configure.ac

@@ -16,7 +16,13 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 # Init the project.
-AC_INIT([SWIFT],[0.4.0])
+AC_INIT([SWIFT],[0.4.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
+swift_config_flags="$*"
+
+#  Need to define this, instead of using fifth argument of AC_INIT, until 2.64.
+AC_DEFINE([PACKAGE_URL],["www.swiftsim.com"], [Package web pages])
+
+AC_COPYRIGHT
 AC_CONFIG_SRCDIR([src/space.c])
 AC_CONFIG_AUX_DIR([.])
 AM_INIT_AUTOMAKE
@@ -182,6 +188,18 @@ if test "$enable_task_debugging" = "yes"; then
    AC_DEFINE([SWIFT_DEBUG_TASKS],1,[Enable task debugging])
 fi
 
+# Check if expensive debugging is on.
+AC_ARG_ENABLE([debugging-checks],
+   [AS_HELP_STRING([--enable-debugging-checks],
+     [Activate expensive consistency checks @<:@yes/no@:>@]
+   )],
+   [enable_debugging_checks="$enableval"],
+   [enable_debugging_checks="no"]
+)
+if test "$enable_debugging_checks" = "yes"; then
+   AC_DEFINE([SWIFT_DEBUG_CHECKS],1,[Enable expensive debugging])
+fi
+
 # Define HAVE_POSIX_MEMALIGN if it works.
 AX_FUNC_POSIX_MEMALIGN
 
@@ -278,7 +296,6 @@ if test "$enable_san" = "yes"; then
    fi
 fi
 
-
 # Autoconf stuff.
 AC_PROG_INSTALL
 AC_PROG_MAKE_SET
@@ -455,8 +472,9 @@ if test "$ac_cv_header_fftw3_h" = "yes"; then
 fi
 AC_SUBST([FFTW_LIBS])
 
-# Check for Intel intrinsics header optionally used by vector.h.
+# Check for Intel and PowerPC intrinsics header optionally used by vector.h.
 AC_CHECK_HEADERS([immintrin.h])
+AC_CHECK_HEADERS([altivec.h])
 
 # Check for timing functions needed by cycle.h.
 AC_HEADER_TIME
@@ -513,6 +531,222 @@ if test "$enable_warn" != "no"; then
     fi
 fi
 
+# Various package configuration options.
+
+# Hydro scheme.
+AC_ARG_WITH([hydro],
+   [AS_HELP_STRING([--with-hydro=<scheme>],
+      [Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo default: gadget2@:>@]
+   )],
+   [with_hydro="$withval"],
+   [with_hydro="gadget2"]
+)
+case "$with_hydro" in
+   gadget2)
+      AC_DEFINE([GADGET2_SPH], [1], [Gadget-2 SPH])
+   ;; 
+   minimal)
+      AC_DEFINE([MINIMAL_SPH], [1], [Minimal SPH])
+   ;; 
+   hopkins)
+      AC_DEFINE([HOPKINS_PE_SPH], [1], [Pressure-Entropy SPH])
+   ;; 
+   default)
+      AC_DEFINE([DEFAULT_SPH], [1], [Default SPH])
+   ;; 
+   gizmo)
+      AC_DEFINE([GIZMO_SPH], [1], [GIZMO SPH])
+   ;; 
+
+   *)
+      AC_MSG_ERROR([Unknown hydrodynamics scheme: $with_hydro])
+   ;;
+esac
+
+# SPH Kernel function
+AC_ARG_WITH([kernel],
+   [AS_HELP_STRING([--with-kernel=<kernel>],
+      [Kernel function to use @<:@cubic-spline, quartic-spline, quintic-spline, wendland-C2, wendland-C4, wendland-C6 default: cubic-spline@:>@]
+   )],
+   [with_kernel="$withval"],
+   [with_kernel="cubic-spline"]
+)
+case "$with_kernel" in
+   cubic-spline)
+      AC_DEFINE([CUBIC_SPLINE_KERNEL], [1], [Cubic spline kernel])
+   ;;    
+   quartic-spline)
+      AC_DEFINE([QUARTIC_SPLINE_KERNEL], [1], [Quartic spline kernel])
+   ;;    
+   quintic-spline)
+      AC_DEFINE([QUINTIC_SPLINE_KERNEL], [1], [Quintic spline kernel])
+   ;;    
+   wendland-C2)
+      AC_DEFINE([WENDLAND_C2_KERNEL], [1], [Wendland-C2 kernel])
+   ;;    
+   wendland-C4)
+      AC_DEFINE([WENDLAND_C4_KERNEL], [1], [Wendland-C4 kernel])
+   ;;    
+   wendland-C6)
+      AC_DEFINE([WENDLAND_C6_KERNEL], [1], [Wendland-C6 kernel])
+   ;;    
+   *)
+      AC_MSG_ERROR([Unknown kernel function: $with_kernel])
+   ;;
+esac
+
+#  Dimensionality of the hydro scheme.
+AC_ARG_WITH([hydro-dimension],
+   [AS_HELP_STRING([--with-hydro-dimension=<dim>],
+      [dimensionality of problem @<:@3/2/1 default: 3@:>@]
+   )],
+   [with_dimension="$withval"],
+   [with_dimension="3"]
+)
+case "$with_dimension" in
+   1)
+      AC_DEFINE([HYDRO_DIMENSION_1D], [1], [1D analysis])
+   ;; 
+   2)
+      AC_DEFINE([HYDRO_DIMENSION_2D], [2], [2D analysis])
+   ;; 
+   3)
+      AC_DEFINE([HYDRO_DIMENSION_3D], [3], [3D analysis])
+   ;; 
+   *)
+      AC_MSG_ERROR([Dimensionality must be 1, 2 or 3])
+   ;;
+esac
+
+#  Equation of state
+AC_ARG_WITH([equation-of-state],
+   [AS_HELP_STRING([--with-equation-of-state=<EoS>],
+      [equation of state @<:@ideal-gas, isothermal-gas default: ideal-gas@:>@]
+   )],
+   [with_eos="$withval"],
+   [with_eos="ideal-gas"]
+)
+case "$with_eos" in
+   ideal-gas)
+      AC_DEFINE([EOS_IDEAL_GAS], [1], [Ideal gas equation of state])
+   ;; 
+   isothermal-gas)
+      AC_DEFINE([EOS_ISOTHERMAL_GAS], [1], [Isothermal gas equation of state])
+   ;; 
+   *)
+      AC_MSG_ERROR([Unknown equation of state: $with_eos])
+   ;;
+esac
+
+#  Adiabatic index
+AC_ARG_WITH([adiabatic-index],
+   [AS_HELP_STRING([--with-adiabatic-index=<gamma>],
+      [adiabatic index @<:@5/3, 7/5, 4/3, 2 default: 5/3@:>@]
+   )],
+   [with_gamma="$withval"],
+   [with_gamma="5/3"]
+)
+case "$with_gamma" in
+   5/3)
+      AC_DEFINE([HYDRO_GAMMA_5_3], [5./3.], [Adiabatic index is 5/3])
+   ;; 
+   7/5)
+      AC_DEFINE([HYDRO_GAMMA_7_5], [7./5.], [Adiabatic index is 7/5])
+   ;; 
+   4/3)
+      AC_DEFINE([HYDRO_GAMMA_4_3], [4./3.], [Adiabatic index is 4/3])
+   ;; 
+   2)
+      AC_DEFINE([HYDRO_GAMMA_2_1], [2.], [Adiabatic index is 2])
+   ;; 
+   *)
+      AC_MSG_ERROR([Unknown adiabatic index: $with_gamma])
+   ;;
+esac
+
+#  Riemann solver
+AC_ARG_WITH([riemann-solver],
+   [AS_HELP_STRING([--with-riemann-solver=<solver>],
+      [riemann solver (gizmo-sph only) @<:@none, exact, trrs, hllc, default: none@:>@]
+   )],
+   [with_riemann="$withval"],
+   [with_riemann="none"]
+)
+case "$with_riemann" in
+   none)
+      AC_DEFINE([RIEMANN_SOLVER_NONE], [1], [No Riemann solver])
+   ;; 
+   exact)
+      AC_DEFINE([RIEMANN_SOLVER_EXACT], [1], [Exact Riemann solver])
+   ;; 
+   trrs)
+      AC_DEFINE([RIEMANN_SOLVER_TRRS], [1], [Two Rarefaction Riemann Solver])
+   ;; 
+   hllc)
+      AC_DEFINE([RIEMANN_SOLVER_HLLC], [1], [Harten-Lax-van Leer-Contact Riemann solver])
+   ;; 
+   *)
+      AC_MSG_ERROR([Unknown Riemann solver: $with_riemann])
+   ;;
+esac
+
+#  Cooling function
+AC_ARG_WITH([cooling],
+   [AS_HELP_STRING([--with-cooling=<function>],
+      [cooling function @<:@none, const-du, const-lambda, grackle default: none@:>@]
+   )],
+   [with_cooling="$withval"],
+   [with_cooling="none"]
+)
+case "$with_cooling" in
+   none)
+      AC_DEFINE([COOLING_NONE], [1], [No cooling function])
+   ;; 
+   const-du)
+      AC_DEFINE([COOLING_CONST_DU], [1], [Const du/dt cooling function])
+   ;; 
+   const-lambda)
+      AC_DEFINE([COOLING_CONST_LAMBDA], [1], [Const Lambda cooling function])
+   ;; 
+   grackle)
+      AC_DEFINE([COOLING_GRACKLE], [1], [Cooling via the grackle library])
+   ;; 
+   *)
+      AC_MSG_ERROR([Unknown cooling function: $with_cooling])
+   ;;
+esac
+
+#  External potential
+AC_ARG_WITH([ext-potential],
+   [AS_HELP_STRING([--with-ext-potential=<pot>],
+      [external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch default: none@:>@]
+   )],
+   [with_potential="$withval"],
+   [with_potential="none"]
+)
+case "$with_potential" in
+   none)
+      AC_DEFINE([EXTERNAL_POTENTIAL_NONE], [1], [No external potential])
+   ;; 
+   point-mass)
+      AC_DEFINE([EXTERNAL_POTENTIAL_POINTMASS], [1], [Point-mass external potential])
+   ;; 
+   isothermal)
+      AC_DEFINE([EXTERNAL_POTENTIAL_ISOTHERMAL], [1], [Isothermal external potential])
+   ;; 
+   softened-isothermal)
+      AC_DEFINE([EXTERNAL_POTENTIAL_SOFTENED_ISOTHERMAL], [1], [Softened isothermal external potential])
+   ;; 
+   disc-patch)
+      AC_DEFINE([EXTERNAL_POTENTIAL_DISC_PATCH], [1], [Disc-patch external potential])
+   ;; 
+   *)
+      AC_MSG_ERROR([Unknown external potential: $with_potential])
+   ;;
+esac
+
+
+
 # Check for git, needed for revision stamps.
 AC_PATH_PROG([GIT_CMD], [git])
 AC_SUBST([GIT_CMD])
@@ -531,6 +765,9 @@ AC_CONFIG_FILES([tests/test27cellsPerturbed.sh], [chmod +x tests/test27cellsPert
 AC_CONFIG_FILES([tests/test125cells.sh], [chmod +x tests/test125cells.sh])
 AC_CONFIG_FILES([tests/testParser.sh], [chmod +x tests/testParser.sh])
 
+# Save the compilation options
+AC_DEFINE_UNQUOTED([SWIFT_CONFIG_FLAGS],["$swift_config_flags"],[Flags passed to configure])
+
 # Report general configuration.
 AC_MSG_RESULT([
    Compiler        : $CC
@@ -545,7 +782,17 @@ AC_MSG_RESULT([
    libNUMA enabled : $have_numa
    Using tcmalloc  : $have_tcmalloc
    CPU profiler    : $have_profiler
-   Task debugging  : $enable_task_debugging
+
+   Hydro scheme       : $with_hydro
+   Dimensionality     : $with_dimension
+   Kernel function    : $with_kernel
+   Equation of state  : $with_eos
+   Adiabatic index    : $with_gamma
+   Riemann solver     : $with_riemann
+   Cooling function   : $with_cooling
+   External potential : $with_potential
+   Task debugging     : $enable_task_debugging
+   Debugging checks   : $enable_debugging_checks
 ])
 
 # Generate output.

examples/CoolingHaloWithSpin/README

-
+Initial Conditions Generation
+-----------------------------
 To make the initial conditions we distribute gas particles randomly in
 a cube with a side length twice that of the virial radius. The density
 profile of the gas is proportional to r^(-2) where r is the distance
@@ -16,14 +17,31 @@ While the system is initially in hydrostatic equilibrium, the cooling
 of the gas and the non-zero angular momentum means that the halo will
 collapse into a spinning disc.
 
+Compilation
+-----------
 To run this example, make such that the code is compiled with either
 the isothermal potential or softened isothermal potential, and
 'const_lambda' cooling, set in src/const.h. In the latter case, a
 (small) value of epsilon needs to be set in cooling.yml.  0.1 kpc
 should work well.
 
+Checking Results
+----------------
 The plotting scripts produce a plot of the density, internal energy
 and radial velocity profile for each
 snapshot. test_energy_conservation.py shows the evolution of energy
 with time. These can be used to check if the example has run properly.
 
+Generating Video
+----------------
+If you want to generate a video of the simulation, the frequency of 
+the snaphots needs to be increased. This can be modified in cooling.yml 
+by changing 'delta_time' to 0.01. 
+
+Once you have the snapshots, 'gadgetviewer' can be used to create a 
+series of snapshot frames. The frames can then be combined together with 
+'ffmpeg' to produce a video. The following command can be used:
+
+ffmpeg -r 20 -i frame_%05d.image.png -c:v ffv1 -qscale:v 0 movie.avi
+
+to produce the video.

examples/CoolingHaloWithSpin/test_energy_conservation.py

@@ -2,8 +2,11 @@ import numpy as np
 import h5py as h5
 import matplotlib.pyplot as plt
 import sys
+import glob
 
-n_snaps = 101
+# Get the total number of snapshots
+file_list = glob.glob("CoolingHalo_*")
+n_snaps = len(file_list)
 
 #some constants
 OMEGA = 0.3 # Cosmological matter fraction at z = 0

src/Makefile.am

@@ -44,8 +44,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h \
     physical_constants.h physical_constants_cgs.h potential.h version.h \
     hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
-    sourceterms_struct.h statistics.h cache.h runner_doiact_vec.h
-
+    sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h
 
 # Common source files
 AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -83,7 +82,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
 		 potential/softened_isothermal/potential.h \
 		 cooling/none/cooling.h cooling/none/cooling_struct.h \
 	         cooling/const_du/cooling.h cooling/const_du/cooling_struct.h \
-                 cooling/const_lambda/cooling.h cooling/const_lambda/cooling_struct.h 
+                 cooling/const_lambda/cooling.h cooling/const_lambda/cooling_struct.h \
+                 memswap.h
 
 
 # Sources and flags for regular library
@@ -108,12 +108,14 @@ version_string.h: version_string.h.in $(AM_SOURCES) $(include_HEADERS) $(noinst_
 	    GIT_BRANCH=`$(GIT_CMD) branch | sed -n 's/^\* \(.*\)/\1/p'`; \
 	    sed -e "s,@PACKAGE_VERSION\@,$(PACKAGE_VERSION)," \
 	        -e "s,@GIT_REVISION\@,$${GIT_REVISION}," \
-	        -e "s|@GIT_BRANCH\@|$${GIT_BRANCH}|" $< > version_string.h; \
+	        -e "s|@GIT_BRANCH\@|$${GIT_BRANCH}|" \
+	        -e "s|@SWIFT_CFLAGS\@|$(CFLAGS)|" $< > version_string.h; \
 	else \
 	    if test ! -f version_string.h; then \
 	        sed -e "s,@PACKAGE_VERSION\@,$(PACKAGE_VERSION)," \
 	            -e "s,@GIT_REVISION\@,unknown," \
-		    -e "s,@GIT_BRANCH\@,unknown," $< > version_string.h; \
+		    -e "s,@GIT_BRANCH\@,unknown," \
+	            -e "s|@SWIFT_CFLAGS\@|$(CFLAGS)|" $< > version_string.h; \
 	    fi; \
 	fi
 

src/active.h

@@ -24,7 +24,6 @@
 
 /* Local includes. */
 #include "cell.h"
-#include "const.h"
 #include "engine.h"
 #include "part.h"
 

src/adiabatic_index.h

@@ -33,7 +33,6 @@
 #include <math.h>
 
 /* Local headers. */
-#include "const.h"
 #include "debug.h"
 #include "error.h"
 #include "inline.h"

src/cell.c

@@ -53,6 +53,8 @@
 #include "gravity.h"
 #include "hydro.h"
 #include "hydro_properties.h"
+#include "memswap.h"
+#include "minmax.h"
 #include "scheduler.h"
 #include "space.h"
 #include "timers.h"
@@ -463,122 +465,73 @@ void cell_gunlocktree(struct cell *c) {
  * @param c The #cell array to be sorted.
  * @param parts_offset Offset of the cell parts array relative to the
  *        space's parts array, i.e. c->parts - s->parts.
+ * @param buff A buffer with at least max(c->count, c->gcount) entries,
+ *        used for sorting indices.
  */
-void cell_split(struct cell *c, ptrdiff_t parts_offset) {
+void cell_split(struct cell *c, ptrdiff_t parts_offset, int *buff) {
 
-  int i, j;
   const int count = c->count, gcount = c->gcount;
   struct part *parts = c->parts;
   struct xpart *xparts = c->xparts;
   struct gpart *gparts = c->gparts;
-  int left[8], right[8];
-  double pivot[3];
-
-  /* Init the pivots. */
-  for (int k = 0; k < 3; k++) pivot[k] = c->loc[k] + c->width[k] / 2;
-
-  /* Split along the x-axis. */
-  i = 0;
-  j = count - 1;
-  while (i <= j) {
-    while (i <= count - 1 && parts[i].x[0] <= pivot[0]) i += 1;
-    while (j >= 0 && parts[j].x[0] > pivot[0]) j -= 1;
-    if (i < j) {
-      struct part temp = parts[i];
-      parts[i] = parts[j];
-      parts[j] = temp;
-      struct xpart xtemp = xparts[i];
-      xparts[i] = xparts[j];
-      xparts[j] = xtemp;
-    }
+  const double pivot[3] = {c->loc[0] + c->width[0] / 2,
+                           c->loc[1] + c->width[1] / 2,
+                           c->loc[2] + c->width[2] / 2};
+  int bucket_count[8] = {0, 0, 0, 0, 0, 0, 0, 0};
+  int bucket_offset[9];
+
+  /* If the buff is NULL, allocate it, and remember to free it. */
+  const int allocate_buffer = (buff == NULL);
+  if (allocate_buffer &&
+      (buff = (int *)malloc(sizeof(int) * max(count, gcount))) == NULL)
+    error("Failed to allocate temporary indices.");
+
+  /* Fill the buffer with the indices. */
+  for (int k = 0; k < count; k++) {
+    const int bid = (parts[k].x[0] > pivot[0]) * 4 +
+                    (parts[k].x[1] > pivot[1]) * 2 + (parts[k].x[2] > pivot[2]);
+    bucket_count[bid]++;
+    buff[k] = bid;
   }
 
-#ifdef SWIFT_DEBUG_CHECKS
-  for (int k = 0; k <= j; k++)
-    if (parts[k].x[0] > pivot[0]) error("cell_split: sorting failed.");
-  for (int k = i; k < count; k++)
-    if (parts[k].x[0] < pivot[0]) error("cell_split: sorting failed.");
-#endif
-
-  left[1] = i;
-  right[1] = count - 1;
-  left[0] = 0;
-  right[0] = j;
-
-  /* Split along the y axis, twice. */
-  for (int k = 1; k >= 0; k--) {
-    i = left[k];
-    j = right[k];
-    while (i <= j) {
-      while (i <= right[k] && parts[i].x[1] <= pivot[1]) i += 1;
-      while (j >= left[k] && parts[j].x[1] > pivot[1]) j -= 1;
-      if (i < j) {
-        struct part temp = parts[i];
-        parts[i] = parts[j];
-        parts[j] = temp;
-        struct xpart xtemp = xparts[i];
-        xparts[i] = xparts[j];
-        xparts[j] = xtemp;
-      }
-    }
-
-#ifdef SWIFT_DEBUG_CHECKS
-    for (int kk = left[k]; kk <= j; kk++)
-      if (parts[kk].x[1] > pivot[1]) {
-        message("ival=[%i,%i], i=%i, j=%i.", left[k], right[k], i, j);
-        error("sorting failed (left).");
-      }
-    for (int kk = i; kk <= right[k]; kk++)
-      if (parts[kk].x[1] < pivot[1]) error("sorting failed (right).");
-#endif
-
-    left[2 * k + 1] = i;
-    right[2 * k + 1] = right[k];
-    left[2 * k] = left[k];
-    right[2 * k] = j;
+  /* Set the buffer offsets. */
+  bucket_offset[0] = 0;
+  for (int k = 1; k <= 8; k++) {
+    bucket_offset[k] = bucket_offset[k - 1] + bucket_count[k - 1];
+    bucket_count[k - 1] = 0;
   }
 
-  /* Split along the z axis, four times. */
-  for (int k = 3; k >= 0; k--) {
-    i = left[k];
-    j = right[k];
-    while (i <= j) {
-      while (i <= right[k] && parts[i].x[2] <= pivot[2]) i += 1;
-      while (j >= left[k] && parts[j].x[2] > pivot[2]) j -= 1;
-      if (i < j) {
-        struct part temp = parts[i];
-        parts[i] = parts[j];
-        parts[j] = temp;
-        struct xpart xtemp = xparts[i];
-        xparts[i] = xparts[j];
-        xparts[j] = xtemp;
+  /* Run through the buckets, and swap particles to their correct spot. */
+  for (int bucket = 0; bucket < 8; bucket++) {
+    for (int k = bucket_offset[bucket] + bucket_count[bucket];
+         k < bucket_offset[bucket + 1]; k++) {
+      int bid = buff[k];
+      if (bid != bucket) {
+        struct part part = parts[k];
+        struct xpart xpart = xparts[k];
+        while (bid != bucket) {
+          int j = bucket_offset[bid] + bucket_count[bid]++;
+          while (buff[j] == bid) {
+            j++;
+            bucket_count[bid]++;
+          }
+          memswap(&parts[j], &part, sizeof(struct part));
+          memswap(&xparts[j], &xpart, sizeof(struct xpart));
+          memswap(&buff[j], &bid, sizeof(int));
+        }
+        parts[k] = part;
+        xparts[k] = xpart;
+        buff[k] = bid;
       }
+      bucket_count[bid]++;
     }
-
-#ifdef SWIFT_DEBUG_CHECKS
-    for (int kk = left[k]; kk <= j; kk++)
-      if (parts[kk].x[2] > pivot[2]) {
-        message("ival=[%i,%i], i=%i, j=%i.", left[k], right[k], i, j);
-        error("sorting failed (left).");
-      }
-    for (int kk = i; kk <= right[k]; kk++)
-      if (parts[kk].x[2] < pivot[2]) {
-        message("ival=[%i,%i], i=%i, j=%i.", left[k], right[k], i, j);
-        error("sorting failed (right).");
-      }
-#endif
-
-    left[2 * k + 1] = i;
-    right[2 * k + 1] = right[k];
-    left[2 * k] = left[k];
-    right[2 * k] = j;
   }
 
   /* Store the counts and offsets. */
   for (int k = 0; k < 8; k++) {
-    c->progeny[k]->count = right[k] - left[k] + 1;
-    c->progeny[k]->parts = &c->parts[left[k]];
-    c->progeny[k]->xparts = &c->xparts[left[k]];
+    c->progeny[k]->count = bucket_count[k];
+    c->progeny[k]->parts = &c->parts[bucket_offset[k]];
+    c->progeny[k]->xparts = &c->xparts[bucket_offset[k]];
   }
 
   /* Re-link the gparts. */
@@ -614,66 +567,51 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset) {
 #endif
 
   /* Now do the same song and dance for the gparts. */
-
-  /* Split along the x-axis. */
-  i = 0;
-  j = gcount - 1;
-  while (i <= j) {
-    while (i <= gcount - 1 && gparts[i].x[0] <= pivot[0]) i += 1;
-    while (j >= 0 && gparts[j].x[0] > pivot[0]) j -= 1;
-    if (i < j) {
-      struct gpart gtemp = gparts[i];
-      gparts[i] = gparts[j];
-      gparts[j] = gtemp;
-    }
+  for (int k = 0; k < 8; k++) bucket_count[k] = 0;
+
+  /* Fill the buffer with the indices. */
+  for (int k = 0; k < gcount; k++) {
+    const int bid = (gparts[k].x[0] > pivot[0]) * 4 +
+                    (gparts[k].x[1] > pivot[1]) * 2 +
+                    (gparts[k].x[2] > pivot[2]);
+    bucket_count[bid]++;
+    buff[k] = bid;
   }
-  left[1] = i;
-  right[1] = gcount - 1;
-  left[0] = 0;
-  right[0] = j;
-
-  /* Split along the y axis, twice. */
-  for (int k = 1; k >= 0; k--) {
-    i = left[k];
-    j = right[k];
-    while (i <= j) {
-      while (i <= right[k] && gparts[i].x[1] <= pivot[1]) i += 1;
-      while (j >= left[k] && gparts[j].x[1] > pivot[1]) j -= 1;
-      if (i < j) {
-        struct gpart gtemp = gparts[i];
-        gparts[i] =