diff --git a/examples/plot_sorted_all.py b/examples/plot_sorted_all.py
new file mode 100644
index 0000000000000000000000000000000000000000..077f5ead53e61c4e85a58094a479b8eebf001b57
--- /dev/null
+++ b/examples/plot_sorted_all.py
@@ -0,0 +1,185 @@
+#!/usr/bin/env python2
+# -*- coding: utf-8 -*-
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import numpy
+# tex stuff
+#rc('font',**{'family':'serif','serif':['Palatino']})
+params = {'axes.labelsize': 16,
+ 'axes.titlesize': 20,
+ 'text.fontsize': 14,
+ 'legend.fontsize': 14,
+ 'xtick.labelsize': 14,
+ 'ytick.labelsize': 14,
+ 'text.usetex': True,
+
+ 'figure.figsize' : (12,9),
+ 'figure.subplot.left' : 0.07, # the left side of the subplots of the figure
+ 'figure.subplot.right' : 0.98 , # the right side of the subplots of the figure
+ 'figure.subplot.bottom' : 0.09 , # the bottom of the subplots of the figure
+ 'figure.subplot.top' : 0.9 , # the top of the subplots of the figure
+ 'figure.subplot.wspace' : 0.26 , # the amount of width reserved for blank space between subplots
+ 'figure.subplot.hspace' : 0.22 , # the amount of height reserved for white space between subplots
+
+ 'lines.markersize' : 2,
+ 'lines.linewidth' : 2,
+
+# 'axes.formatter.limits' : (-2, 1),
+
+ 'text.latex.unicode': True
+ }
+rcParams.update(params)
+rc('font', family='serif')
+import sys
+import os
+
+print "Plotting..."
+
+
+
+# Read Quickshed accelerations
+data=loadtxt( "interaction_dump_all.dat" )
+id = data[:,0]
+
+mass = data[:,1]
+pos_x = data[:,2]
+pos_y = data[:,3]
+pos_z = data[:,4]
+
+accx_u=data[:,5]
+accy_u=data[:,6]
+accz_u=data[:,7]
+
+accx_s=data[:,8]
+accy_s=data[:,9]
+accz_s=data[:,10]
+
+
+# Build error ------------------------------------------------
+
+errx_s = (accx_s - accx_u )/sqrt(accx_u**2 + accy_u**2 + accz_u**2)
+erry_s = (accy_s - accy_u )/sqrt(accy_u**2 + accy_u**2 + accz_u**2)
+errz_s = (accz_s - accz_u )/sqrt(accx_u**2 + accy_u**2 + accz_u**2)
+
+e_errx_s = errx_s#[abs(errx_s) > 0.001]
+e_erry_s = erry_s#[abs(erry_s) > 0.001]
+e_errz_s = errz_s#[abs(errz_s) > 0.001]
+
+# Statistics
+meanx_s = mean(errx_s[abs(errx_s) < 0.1])
+stdx_s = std(errx_s[abs(errx_s) < 0.1])
+meany_s = mean(erry_s[abs(erry_s) < 0.1])
+stdy_s = std(erry_s[abs(erry_s) < 0.1])
+meanz_s = mean(errz_s[abs(errz_s) < 0.1])
+stdz_s = std(errz_s[abs(errz_s) < 0.1])
+
+
+
+# Plot -------------------------------------------------------
+figure(frameon=True)
+
+subplot(311, title="Acceleration along X")
+#plot(id[abs(errx_s) > 0.001], e_errx_s , 'ro')
+plot(pos_x, e_errx_s , 'ro')
+#text( 0., 0.1, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
+ylim(-0.2, 0.2)
+grid()
+
+subplot(312, title="Acceleration along Y")
+#plot(id[abs(erry_s) > 0.001], e_erry_s , 'ro')
+plot(pos_y, e_erry_s , 'ro')
+#text( 0., 0.1, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
+ylim(-0.2, 0.2)
+grid()
+
+subplot(313, title="Acceleration along Z")
+#plot(id[abs(errz_s) > 0.001], e_errz_s , 'ro', label="Sorted")
+plot(pos_z, e_errz_s , 'ro', label="Sorted")
+#text( 0., 0.1, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
+
+legend(loc="upper right")
+
+ylim(-0.2, 0.2)
+#xlim(0, size(id)-1)
+grid()
+
+savefig("accelerations_relative_all.png")
+close()
+
+
+# Plot -------------------------------------------------------
+figure(frameon=True)
+
+subplot(311, title="Acceleration along X")
+#plot(id[abs(errx_s) > 0.001], e_errx_s , 'ro')
+plot(pos_x, accx_u - accx_s, 'bx')
+#plot(pos_x, accx_s , 'ro')
+#text( 0., 0.1, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
+ylim(-70, 70)
+grid()
+
+subplot(312, title="Acceleration along Y")
+#plot(id[abs(erry_s) > 0.001], e_erry_s , 'ro')
+plot(pos_y, accy_u - accy_s, 'bx')
+#plot(pos_y, accy_s , 'ro')
+#text( 0., 0.1, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
+ylim(-70, 70)
+grid()
+
+subplot(313, title="Acceleration along Z")
+#plot(id[abs(errz_s) > 0.001], e_errz_s , 'ro', label="Sorted")
+plot(pos_z, accz_u - accz_s , 'bx', label="Unsorted")
+#plot(pos_z, accz_s , 'ro', label="Sorted")
+
+#text( 0., 0.1, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
+
+legend(loc="upper right")
+
+ylim(-70, 70)
+grid()
+
+savefig("accelerations_absolute_all.png")
+close()
+
+
+
+# Plot -------------------------------------------------------
+bins = linspace(-3, 3, 10000)
+
+e_errx_s = errx_s[(abs(errx_s) >= 0.000) & (abs(errx_s) < 1.)]
+e_erry_s = erry_s[(abs(erry_s) >= 0.000) & (abs(erry_s) < 1.)]
+e_errz_s = errz_s[(abs(errz_s) >= 0.000) & (abs(errz_s) < 1.)]
+
+
+
+figure(frameon=True)
+subplot(311, title="Acceleration along X")
+hist(e_errx_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r', label="Sorted")
+legend(loc="upper right")
+xlim(-0.12, 0.12)
+
+subplot(312, title="Acceleration along Y")
+hist(e_erry_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
+xlim(-0.12, 0.12)
+
+subplot(313, title="Acceleration along Z")
+hist(e_errz_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r')
+xlim(-0.12, 0.12)
+
+savefig("histogram_all.png")
+close()
+
+
+
+
+print "E_error for sorted: ( x= %4.3e"%std(e_errx_s), "y= %4.3e"%std(e_erry_s), "z= %4.3e"%std(e_errz_s), ") Combined YZ = %4.3e"%(( std(e_erry_s) + std(e_errz_s) )/2.)
+
+ #print std(e_errx_s), std(e_erry_s), std(e_errz_s)
+
+
+ #print "Min for sorted: ( x= %4.3e"%min(errx_s), "y= %4.3e"%min(erry_s), "z= %4.3e"%min(errz_s), ") Combined YZ = %4.3e"%(min( min(erry_s), min(errz_s) ))
+ #print "Max for sorted: ( x= %4.3e"%max(errx_s), "y= %4.3e"%max(erry_s), "z= %4.3e"%max(errz_s), ") Combined YZ = %4.3e"%(max( max(erry_s), max(errz_s) ))
+
+
+
diff --git a/examples/test_bh_sorted.c b/examples/test_bh_sorted.c
index 7f74823e13807f2b7503846a029dbf192ceeff74..f3560c701aa4351e9e6e0ebe9df7c2be37ba95b2 100644
--- a/examples/test_bh_sorted.c
+++ b/examples/test_bh_sorted.c
@@ -40,7 +40,7 @@
#define cell_maxparts 100
#define task_limit 1e8
#define const_G 1 // 6.6738e-8
-#define dist_min 0.5 /* Used fpr legacy walk only */
+#define dist_min 0.5 /* Used for legacy walk only */
#define dist_cutoff_ratio 1.5
#define iact_pair_direct iact_pair_direct_sorted
@@ -482,11 +482,11 @@ void get_axis(struct cell **ci, struct cell **cj, struct index **ind_i,
orth1[3] = -d1;
orth2[3] = -d2;
- message("dx=[%.3e,%.3e,%.3e], axis=[%.3e,%.3e,%.3e]",
- dx[0], dx[1], dx[2], axis[0], axis[1], axis[2]);
- message("N_planes= %d", *num_orth_planes);
- message("o1=[%.3e,%.3e,%.3e] %.3e", orth1[0], orth1[1], orth1[2], orth1[3]);
- message("o2=[%.3e,%.3e,%.3e] %.3e", orth2[0], orth2[1], orth2[2], orth2[3]);
+ /* message("dx=[%.3e,%.3e,%.3e], axis=[%.3e,%.3e,%.3e]", */
+ /* dx[0], dx[1], dx[2], axis[0], axis[1], axis[2]); */
+ /* message("N_planes= %d", *num_orth_planes); */
+ /* message("o1=[%.3e,%.3e,%.3e] %.3e", orth1[0], orth1[1], orth1[2], orth1[3]); */
+ /* message("o2=[%.3e,%.3e,%.3e] %.3e", orth2[0], orth2[1], orth2[2], orth2[3]); */
/* Make sure the sorts are ok. */
/* for (int k = 1; k < (*ci)->count; k++)
@@ -2058,6 +2058,39 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
free(parts);
}
+
+/* All 26 configurations for the cell tests */
+const float cell_shift[27 * 3] = {
+ 1.0, 1.0, 1.0, /* 0 */
+ 1.0, 1.0, 0.0, /* 1 */
+ 1.0, 1.0, -1.0, /* 2 */
+ 1.0, 0.0, 1.0, /* 3 */
+ 1.0, 0.0, 0.0, /* 4 */
+ 1.0, 0.0, -1.0, /* 5 */
+ 1.0, -1.0, 1.0, /* 6 */
+ 1.0, -1.0, 0.0, /* 7 */
+ 1.0, -1.0, -1.0, /* 8 */
+ 0.0, 1.0, 1.0, /* 9 */
+ 0.0, 1.0, 0.0, /* 10 */
+ 0.0, 1.0, -1.0, /* 11 */
+ 0.0, 0.0, 1.0, /* 12 */
+ -1.0, -1.0, -1.0, /* 13 */
+ -1.0, -1.0, 0.0, /* 14 */
+ -1.0, -1.0, 1.0, /* 15 */
+ -1.0, 0.0, -1.0, /* 16 */
+ -1.0, 0.0, 0.0, /* 17 */
+ -1.0, 0.0, 1.0, /* 18 */
+ -1.0, 1.0, -1.0, /* 19 */
+ -1.0, 1.0, 0.0, /* 20 */
+ -1.0, 1.0, 1.0, /* 21 */
+ 0.0, -1.0, -1.0, /* 22 */
+ 0.0, -1.0, 0.0, /* 23 */
+ 0.0, -1.0, 1.0, /* 24 */
+ 0.0, 0.0, -1.0, /* 25 */
+ 0.0, 0.0, 0.0 /* 26 */ /* <-- The cell itself */
+};
+
+
/**
* @brief Creates two neighbouring cells wiht N_parts per celland makes them
* interact using both the
@@ -2072,36 +2105,6 @@ void test_direct_neighbour(int N_parts, int orientation) {
struct part *parts;
struct cell left, right;
- /* All 13 configurations */
- const float cell_shift[26 * 3] = {
- 1.0, 1.0, 1.0, /* 0 */
- 1.0, 1.0, 0.0, /* 1 */
- 1.0, 1.0, -1.0, /* 2 */
- 1.0, 0.0, 1.0, /* 3 */
- 1.0, 0.0, 0.0, /* 4 */
- 1.0, 0.0, -1.0, /* 5 */
- 1.0, -1.0, 1.0, /* 6 */
- 1.0, -1.0, 0.0, /* 7 */
- 1.0, -1.0, -1.0, /* 8 */
- 0.0, 1.0, 1.0, /* 9 */
- 0.0, 1.0, 0.0, /* 10 */
- 0.0, 1.0, -1.0, /* 11 */
- 0.0, 0.0, 1.0, /* 12 */
- -1.0, -1.0, -1.0, /* 13 */
- -1.0, -1.0, 0.0, /* 14 */
- -1.0, -1.0, 1.0, /* 15 */
- -1.0, 0.0, -1.0, /* 16 */
- -1.0, 0.0, 0.0, /* 17 */
- -1.0, 0.0, 1.0, /* 18 */
- -1.0, 1.0, -1.0, /* 19 */
- -1.0, 1.0, 0.0, /* 20 */
- -1.0, 1.0, 1.0, /* 21 */
- 0.0, -1.0, -1.0, /* 22 */
- 0.0, -1.0, 0.0, /* 23 */
- 0.0, -1.0, 1.0, /* 24 */
- 0.0, 0.0, -1.0 /* 25 */
- };
-
if ( orientation >= 26 )
error( "Wrong orientation !" );
@@ -2254,6 +2257,146 @@ void test_direct_neighbour(int N_parts, int orientation) {
free(parts);
}
+
+
+
+
+
+/**
+ * @brief Creates 27 cells and makes the particles in the central one
+ * interact with the other cells
+ * using both the sorted and unsorted interactions
+ * Outputs then the two sets of accelerations for accuracy tests.
+ * @param N_parts Number of particles in each cell
+ */
+void test_all_direct_neighbours(int N_parts) {
+
+ int k, j;
+ struct part *parts;
+ struct cell cells[27];
+
+
+ /* Init and fill the particle array. */
+ if ((parts = (struct part *)malloc(sizeof(struct part) * N_parts * 27)) ==
+ NULL)
+ error("Failed to allocate particle buffer.");
+
+
+ for (j = 0; j < 27 ; ++j) {
+
+ float shift[3];
+ for ( k = 0 ; k < 3 ; ++k )
+ shift[k] = cell_shift[ 3 * j + k ];
+
+ /* Create random set of particles in all cells */
+ for (k = 0; k < N_parts; k++) {
+ parts[j*N_parts + k].id = j*N_parts + k;
+
+ parts[j*N_parts + k].x[0] = ((double)rand()) / RAND_MAX + shift[0];
+ parts[j*N_parts + k].x[1] = ((double)rand()) / RAND_MAX + shift[1];
+ parts[j*N_parts + k].x[2] = ((double)rand()) / RAND_MAX + shift[2];
+
+ parts[j*N_parts + k].mass = ((double)rand()) / RAND_MAX;
+ parts[j*N_parts + k].a[0] = 0.0;
+ parts[j*N_parts + k].a[1] = 0.0;
+ parts[j*N_parts + k].a[2] = 0.0;
+ }
+
+
+ /* Get the cell geometry right */
+ cells[j].loc[0] = shift[0];
+ cells[j].loc[1] = shift[1];
+ cells[j].loc[2] = shift[2];
+ cells[j].h = 1.;
+
+ /* Put the particles in the cell */
+ cells[j].parts = parts + N_parts * j;
+ cells[j].count = N_parts;
+
+ /* Get the linked list right (functions should not recurse but hey...) */
+ cells[j].firstchild = NULL;
+ cells[j].sibling = NULL;
+ cells[j].split = 0;
+
+ /* Get the multipoles right (should also be useless) */
+ comp_com(&cells[j]);
+
+ /* Nothing has been sorted yet */
+ cells[j].indices = NULL;
+ cells[j].sorted = 0;
+ }
+
+
+#ifdef COUNTERS
+ count_direct_unsorted = 0;
+ count_direct_sorted_pp = 0;
+ count_direct_sorted_pm_i = 0;
+ count_direct_sorted_pm_j = 0;
+#endif
+
+ /* Do the interactions without sorting */
+ for (j = 0; j < 26 ; ++j)
+ iact_pair_direct_unsorted(&cells[26], &cells[j]);
+
+ message("Unsorted interactions done ");
+
+ /* Store accelerations */
+ for (k = 0; k < 27*N_parts; k++) {
+ parts[k].a_exact[0] = parts[k].a[0];
+ parts[k].a_exact[1] = parts[k].a[1];
+ parts[k].a_exact[2] = parts[k].a[2];
+ parts[k].a[0] = 0.0;
+ parts[k].a[1] = 0.0;
+ parts[k].a[2] = 0.0;
+ }
+
+ /* Do the interactions with sorting */
+ for (j = 0; j < 26 ; ++j)
+ iact_pair_direct_sorted(&cells[26], &cells[j]);
+
+ message("Sorted interactions done ");
+
+
+ /* Now, output everything */
+ char fileName[100];
+ sprintf(fileName, "interaction_dump_all.dat");
+ message("Writing file '%s'", fileName);
+ FILE *file = fopen(fileName, "w");
+ fprintf(file,
+ "# ID m x y z a_u.x a_u.y a_u.z a_s.x a_s.y a_s.z\n");
+ for (k = 0; k < N_parts; k++) {
+ fprintf(file, "%d %e %e %e %e %e %e %e %e %e %e\n", parts[26*N_parts + k].id,
+ parts[26*N_parts + k].mass, parts[26*N_parts + k].x[0], parts[26*N_parts + k].x[1],
+ parts[26*N_parts + k].x[2], parts[26*N_parts + k].a_exact[0], parts[26*N_parts + k].a_exact[1],
+ parts[26*N_parts + k].a_exact[2], parts[26*N_parts + k].a[0], parts[26*N_parts + k].a[1], parts[26*N_parts + k].a[2]);
+ }
+ fclose(file);
+
+
+
+#ifdef COUNTERS
+ message("Unsorted intereactions: %d", count_direct_unsorted);
+ message("Sorted intereactions PP: %d", count_direct_sorted_pp);
+ message("Sorted intereactions PM (part i): %d", count_direct_sorted_pm_i);
+ message("Sorted intereactions PM (part j): %d", count_direct_sorted_pm_j);
+ message("Sorted intereactions total: %d",
+ count_direct_sorted_pm_j + count_direct_sorted_pm_i +
+ count_direct_sorted_pp);
+// message( "%f %d %d %d %d\n", dist_cutoff_ratio, count_direct_unsorted,
+// count_direct_sorted_pp, count_direct_sorted_pm_i, count_direct_sorted_pm_j );
+#endif
+
+ /* Clean up */
+ free(parts);
+
+}
+
+
+
+
+
+
+
/**
* @brief Main function.
*/
@@ -2328,6 +2471,14 @@ int main(int argc, char *argv[]) {
for ( k = 0 ; k < 26 ; ++k )
test_direct_neighbour(N_parts, k);
+ k++;
+ /* Dump arguments */
+ message("Interacting one cell with %d particles with its 27 neighbours",
+ N_parts);
+
+ test_all_direct_neighbours(N_parts);
+
+
} else {
/* Dump arguments. */