diff --git a/examples/plot_sorted.py b/examples/plot_sorted.py new file mode 100644 index 0000000000000000000000000000000000000000..eda2093cb235ecac5214e52be54c9de8d2c82073 --- /dev/null +++ b/examples/plot_sorted.py @@ -0,0 +1,122 @@ +#!/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' : 6, + '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.dat") +id = data[:,0] +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 )/abs(accx_u) +erry_s = (accy_s - accy_u )/abs(accy_u) +errz_s = (accz_s - accz_u )/abs(accz_u) + +# 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, errx_s , 'rx') +#text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanx_bh, stdx_bh, meanx_new, stdx_new), backgroundcolor="w", va="top", ha="right" ) +ylim(-0.2, 0.2) +xlim(0, size(id)-1) +grid() + +subplot(312, title="Acceleration along Y") +plot(id, erry_s , 'rx') +#text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meany_bh, stdy_bh, meany_new, stdy_new), backgroundcolor="w", va="top", ha="right" ) +ylim(-0.2, 0.2) +xlim(0, size(id)-1) + +grid() + +subplot(313, title="Acceleration along Z") +plot(id, errz_s , 'rx', label="QuickShed") +#text(id[-1], 0.18, "B-H: $%5.3f\\pm%5.3f$\n QuickShed: $%5.3f\\pm%5.3f$"%(meanz_bh, stdz_bh, meanz_new, stdz_new), backgroundcolor="w", va="top", ha="right" ) +legend(loc="upper right") + +ylim(-0.2, 0.2) +xlim(0, size(id)-1) +grid() + +savefig("accelerations.png") + + + + +# Plot ------------------------------------------------------- +# bins = linspace(-3, 3, 10000) + + +# figure(frameon=True) +# subplot(311, title="Acceleration along X")#, yscale='log') +# hist(errx_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r', label="Legacy") +# legend(loc="upper right") +# xlim(-0.03, 0.03) + +# subplot(312, title="Acceleration along Y") +# hist(erry_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r') +# xlim(-0.03, 0.03) + +# subplot(313, title="Acceleration along Z") +# hist(errz_s, bins=bins, normed=1, histtype='step', rwidth=0.01, color='r') +# xlim(-0.03, 0.03) + +# savefig("histogram.png") + + + diff --git a/examples/test_bh_sorted.c b/examples/test_bh_sorted.c index 38080b62c4fbb01c742ec3e7fbead6396c0df5f7..5516e1cb07545397fb7938ae8238a6432f2b4a9d 100644 --- a/examples/test_bh_sorted.c +++ b/examples/test_bh_sorted.c @@ -85,14 +85,14 @@ struct cell { /* We keep both CoMs and masses to make sure the comp_com calculation is * correct (use an anonymous union to keep variable names compact). */ - // union { + union { /* Information for the legacy walk */ struct multipole legacy; /* Information for the QuickShed walk */ struct multipole new; - // }; + }; int res, com_tid; struct index *indices; @@ -1938,8 +1938,116 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) { /* Clean up. */ qsched_free(&s); + free(parts); +} + + +/** + * @brief Creates two neighbouring cells wiht N_parts per celland makes them interact using both the + * sorted and unsortde interactions. Outputs then the tow sets of accelerations for accuracy tests. + */ +void test_direct_neighbour( int N_parts ) { + + int k; + struct part *parts; + struct cell left, right; + + /* Init and fill the particle array. */ + if ((parts = (struct part *)malloc(sizeof(struct part) * N_parts * 2)) == NULL) + error("Failed to allocate particle buffer."); + + /* Create random set of particles in both cells */ + for (k = 0; k < 2*N_parts; k++) { + parts[k].id = k; + parts[k].x[0] = ((double)rand()) / RAND_MAX + (k >= N_parts ? 1 : 0. ); + parts[k].x[1] = ((double)rand()) / RAND_MAX; + parts[k].x[2] = ((double)rand()) / RAND_MAX; + parts[k].mass = ((double)rand()) / RAND_MAX; + parts[k].a[0] = 0.0; + parts[k].a[1] = 0.0; + parts[k].a[2] = 0.0; + } + + /* Get the cell geometry right */ + left.loc[0] = 0.; left.loc[1] = 0.; left.loc[2] = 0.; + left.h = 1.; + right.loc[0] = 1.; right.loc[1] = 0.; right.loc[2] = 0.; + right.h = 1.; + + /* Put the particles in the cell */ + left.parts = parts; + left.count = N_parts; + right.parts = parts + N_parts; + right.count = N_parts; + + /* Get the linked list right (functions should not recurse but hey...) */ + left.firstchild = NULL; left.sibling = NULL; + left.split = 0; + right.firstchild = NULL; right.sibling = NULL; + right.split = 0; + + /* Get the multipoles right (should also be useless) */ + comp_com( &left ); + comp_com( &right ); + + /* Nothing has been sorted yet */ + left.indices = NULL; + left.sorted = 0; + right.indices = NULL; + right.sorted = 0; + + + /* Do the interactions without sorting */ + iact_pair_direct_unsorted( &left, &right ); + + message( "Unsorted interactions done " ); + + + /* Store accelerations */ + for (k = 0; k < 2*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 */ + iact_pair_direct_sorted( &left, &right ); + + message( "Sorted interactions done " ); + + + + /* Sort the particles along axis to simply interpretation */ + int compParts(const void* c1, const void* c2) { + if ( ((struct part*)c1)->x[0] < ((struct part*)c2)->x[0] ) return -1; + else if ( ((struct part*)c1)->x[0] == ((struct part*)c2)->x[0] ) return 0; + else return 1; + } + qsort( parts, 2*N_parts, sizeof(struct part), compParts); + + + /* Now, output everything */ + message( "Writing file 'interaction_dump.dat'" ); + FILE* file = fopen("interaction_dump.dat", "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 < 2*N_parts; k++) { + fprintf(file, "%d %e %e %e %e %e %e %e %e %e %e\n", parts[k].id, parts[k].mass, + parts[k].x[0], parts[k].x[1], parts[k].x[2], + parts[k].a_exact[0], parts[k].a_exact[1], parts[k].a_exact[2], + parts[k].a[0], parts[k].a[1], parts[k].a[2]); + } + fclose( file ); + + /* Clean up */ + free( parts ); } + /** * @brief Main function. */ @@ -1949,6 +2057,7 @@ int main(int argc, char *argv[]) { int c, nr_threads; int N = 1000, runs = 1; char fileName[100] = {0}; + int N_parts = 0; /* Die on FP-exceptions. */ feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); @@ -1960,7 +2069,7 @@ int main(int argc, char *argv[]) { } /* Parse the options */ - while ((c = getopt(argc, argv, "n:r:t:f:")) != -1) switch (c) { + while ((c = getopt(argc, argv, "n:r:t:f:c:")) != -1) switch (c) { case 'n': if (sscanf(optarg, "%d", &N) != 1) error("Error parsing number of particles."); @@ -1978,14 +2087,19 @@ int main(int argc, char *argv[]) { if (sscanf(optarg, "%s", &fileName[0]) != 1) error("Error parsing file name."); break; + case 'c': + if (sscanf(optarg, "%d", &N_parts) != 1) + error("Error parsing number of particles in neighbouring cells test."); + break; case '?': fprintf(stderr, - "Usage: %s [-t nr_threads] [-n N] [-r runs] [-f file]\n", + "Usage: %s [-t nr_threads] [-n N] [-r runs] [-f file] [-c Nparts]\n", argv[0]); fprintf(stderr, "Solves the N-body problem using a Barnes-Hut\n" "tree code with N random particles read from a file in " - "[0,1]^3 using\n" - "nr_threads threads.\n"); + "[0,1]^3 using" + "nr_threads threads.\n" + "A test of the neighbouring cells interaction with Nparts per cell is also run.\n" ); exit(EXIT_FAILURE); } @@ -1995,19 +2109,33 @@ int main(int argc, char *argv[]) { /* Part information */ printf("Size of part: %zu bytes.\n", sizeof(struct part)); - /* Dump arguments. */ - if (fileName[0] == 0) { - message("Computing the N-body problem over %i random particles using %i " - "threads (%i runs).", - N, nr_threads, runs); - } else { - message("Computing the N-body problem over %i particles read from '%s' " - "using %i threads (%i runs).", - N, fileName, nr_threads, runs); + /* Run the neighbour direct integration test */ + if ( N_parts > 0 ) { + + /* Dump arguments */ + message("Interacting 2 neighbouring cells with %d particles per cell", N_parts); + + /* Run the test */ + test_direct_neighbour( N_parts ); + } + else { - /* Run the test. */ - test_bh(N, nr_threads, runs, fileName); + /* Dump arguments. */ + if (fileName[0] == 0) { + message("Computing the N-body problem over %i random particles using %i " + "threads (%i runs).", + N, nr_threads, runs); + } else { + message("Computing the N-body problem over %i particles read from '%s' " + "using %i threads (%i runs).", + N, fileName, nr_threads, runs); + } + + /* Run the BH test. */ + test_bh(N, nr_threads, runs, fileName); + + } return 0;