diff --git a/examples/plot_sorted.py b/examples/plot_sorted.py
index e76f338b56d3897d547599f0a23c4fc2fa6e9767..9e07de4edc6d0afad13067742e1a929ed21c8b05 100644
--- a/examples/plot_sorted.py
+++ b/examples/plot_sorted.py
@@ -34,9 +34,12 @@ rc('font', family='serif')
import sys
import os
from scipy import stats
+from scipy import optimize
-dist_cutoff_ratio=1.5
-dist_max = 0.5
+dist_cutoff_ratio=1.2
+ortho_max_fit = 0.2
+ortho_min_fit = 0.
+ortho_max = 2.5
plot_range = 0.01
print "Plotting..."
@@ -73,11 +76,11 @@ axis = [
# Compute correctness limit
limit_exact = ( dist_cutoff_ratio - 1. )
-
+print limit_exact
#names = ["side", "edge", "corner"]
#for orientation in range( 26 ):
-for jjj in range(3):
+for jjj in range(0,3):
if jjj == 0:
orientation = 0
if jjj == 1:
@@ -96,7 +99,8 @@ for jjj in range(3):
y = data[:,4]
z = data[:,5]
- dist = data[:,12]
+ ortho = data[:,12]
+ dist = data[:,13]
accx_u=data[:,6]
accy_u=data[:,7]
@@ -105,16 +109,107 @@ for jjj in range(3):
accx_s=data[:,9]
accy_s=data[:,10]
accz_s=data[:,11]
+
+ #figure()
+ #plot(x, dist,'bx')
+ #savefig("test.png")
+
+ # # Now, prepare the correction ----------------------------------
+ def parabola(x,a,b,c):
+ return a*x**2 + b*x + c
+ def cubic(x,a,b,c,d):
+ return a*x**3 + b*x**2 + c*x + d
+ def quartic(x,a,b,c,d,e):
+ return a*x**4 + b*x**3 + c*x**2 + d*x + e
+
+ def quintic(x,a,b,c,d,e,f):
+ return a*x**5 + b*x**4 + c*x**3 + d*x**2 + e*x + f
+
+
+ cell_edge = sqrt(sum(abs(array(axis[orientation*3:orientation*3+3]))))
+ print cell_edge
+
+ figure()
+
+ xx = dist[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
+ sim = accx_s[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
+ exact = accx_u[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
+
+ sim = sim[xx > cell_edge]
+ exact = exact[xx > cell_edge]
+ xx = xx[xx > cell_edge]
+
+ sim = sim[abs(cell_edge-xx) > limit_exact]
+ exact = exact[abs(cell_edge-xx) > limit_exact]
+ xx = xx[abs(cell_edge-xx) > limit_exact]
+
+ delta = sim / exact
+ parabola_param1,pcov = optimize.curve_fit(parabola, xx, delta, p0=[1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
+ print parabola_param1
+ cubic_param1,pcov = optimize.curve_fit(cubic, xx, delta, p0=[1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
+ print cubic_param1
+ quartic_param1,pcov = optimize.curve_fit(quartic, xx, delta, p0=[1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
+ print quartic_param1
+ quintic_param1,pcov = optimize.curve_fit(quintic, xx, delta, p0=[1,1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
+ print quintic_param1
+
+ plot(xx, delta, 'bx')
+ plot([0,2*cell_edge], [1,1], 'r')
+ xxx = linspace(cell_edge+limit_exact, 2*cell_edge)
+ plot(xxx, parabola(xxx, parabola_param1[0], parabola_param1[1], parabola_param1[2]), 'k-')
+ plot(xxx, cubic(xxx, cubic_param1[0], cubic_param1[1], cubic_param1[2], cubic_param1[3]), 'k--')
+ plot(xxx, quartic(xxx, quartic_param1[0], quartic_param1[1], quartic_param1[2], quartic_param1[3], quartic_param1[4]), 'k:')
+ plot(xxx, quintic(xxx, quintic_param1[0], quintic_param1[1], quintic_param1[2], quintic_param1[3], quintic_param1[4], quintic_param1[5]), 'k-.')
+
+ xx = dist[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
+ sim = accx_s[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
+ exact = accx_u[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
+
+ sim = sim[xx < cell_edge]
+ exact = exact[xx < cell_edge]
+ xx = xx[xx < cell_edge]
+
+ sim = sim[abs(cell_edge-xx) > limit_exact]
+ exact = exact[abs(cell_edge-xx) > limit_exact]
+ xx = xx[abs(cell_edge-xx) > limit_exact]
+
+ delta = sim / exact
+ parabola_param2,pcov = optimize.curve_fit(parabola, xx, delta, p0=[1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
+ print parabola_param2
+ cubic_param2,pcov = optimize.curve_fit(cubic, xx, delta, p0=[1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
+ print cubic_param2
+ quartic_param2,pcov = optimize.curve_fit(quartic, xx, delta, p0=[1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
+ print quartic_param2
+ quintic_param2,pcov = optimize.curve_fit(quintic, xx, delta, p0=[1,1,1,1,1,1], sigma=None, maxfev=100000, xtol=1e-13, ftol=1e-13)
+ print quintic_param2
+
+ plot(xx, delta, 'bx')
+ xxx = linspace(0., cell_edge-limit_exact)
+ plot(xxx, parabola(xxx, parabola_param2[0], parabola_param2[1], parabola_param2[2]), 'k-')
+ plot(xxx, cubic(xxx, cubic_param2[0], cubic_param2[1], cubic_param2[2], cubic_param2[3]), 'k--')
+ plot(xxx, quartic(xxx, quartic_param2[0], quartic_param2[1], quartic_param2[2], quartic_param2[3], quartic_param2[4]), 'k:')
+ plot(xxx, quintic(xxx, quintic_param2[0], quintic_param2[1], quintic_param2[2], quintic_param2[3], quintic_param2[4], quintic_param2[5]), 'k-.')
+
+ savefig("fit_%d.png"%orientation)
+
+ # Apply the correction ---------------------------------------
+ for i in range(size(x)):
+ if cell_edge-dist[i] > limit_exact:
+ accx_s[i] /= quintic(dist[i], quintic_param2[0], quintic_param2[1], quintic_param2[2], quintic_param2[3], quintic_param2[4], quintic_param2[5])
+ if dist[i]-cell_edge > limit_exact:
+ accx_s[i] /= quintic(dist[i], quintic_param1[0], quintic_param1[1], quintic_param1[2], quintic_param1[3], quintic_param1[4], quintic_param1[5])
+
+
# Build error ------------------------------------------------
- inv_acc_tot = 1.0 / sqrt(accx_u**2 + accy_u**2 + accz_u**2)
- errx_s = (accx_s - accx_u ) * inv_acc_tot
- erry_s = (accy_s - accy_u ) * inv_acc_tot
- errz_s = (accz_s - accz_u ) * inv_acc_tot
+ inv_acc_tot = 1.0 / sqrt(accx_u**2 + accy_u**2 + accz_u**2)
+ errx_s = (accx_s - accx_u ) * inv_acc_tot # / abs(accx_u)
+ erry_s = (accy_s - accy_u ) * inv_acc_tot #/ abs(accy_u)
+ errz_s = (accz_s - accz_u ) * inv_acc_tot #/ abs(accz_u)
e_errx_s = errx_s#[abs(errx_s) > 0.001]
e_erry_s = erry_s#[abs(erry_s) > 0.001]
@@ -129,10 +224,10 @@ for jjj in range(3):
stdz_s = std(errz_s[abs(errz_s) < 0.1])
-# sample_pos = pos[dist<0.2]
-# sample_x = e_errx_s[dist<0.2]
-# sample_y = e_erry_s[dist<0.2]
-# sample_z = e_errz_s[dist<0.2]
+# sample_pos = pos[ortho<0.2]
+# sample_x = e_errx_s[ortho<0.2]
+# sample_y = e_erry_s[ortho<0.2]
+# sample_z = e_errz_s[ortho<0.2]
# numBins = 100
@@ -157,7 +252,8 @@ for jjj in range(3):
subplot(311, title="Acceleration along X")
- scatter(pos[dist<dist_max], e_errx_s[dist<dist_max] ,c=dist[dist<dist_max], marker='o', s=1, linewidth=0, cmap='jet')
+ scatter(pos[ortho<ortho_max], e_errx_s[ortho<ortho_max] ,c=ortho[ortho<ortho_max], marker='o', s=1, linewidth=0, cmap='jet')
+ #scatter(pos, e_errx_s ,c=ortho, marker='o', s=1, linewidth=0, cmap='jet')
plot([-limit_exact, -limit_exact], [-plot_range, plot_range], 'k--')
plot([limit_exact, limit_exact], [-plot_range, plot_range], 'k--')
text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
@@ -166,7 +262,7 @@ for jjj in range(3):
grid()
subplot(312, title="Acceleration along Y")
- scatter(pos[dist<dist_max], e_erry_s[dist<dist_max] , c=dist[dist<dist_max], marker='o', s=1, linewidth=0, cmap='jet')
+ scatter(pos[ortho<ortho_max], e_erry_s[ortho<ortho_max] , c=ortho[ortho<ortho_max], marker='o', s=1, linewidth=0, cmap='jet')
plot([-limit_exact, -limit_exact], [-plot_range, plot_range], 'k--')
plot([limit_exact, limit_exact], [-plot_range, plot_range], 'k--')
text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
@@ -175,7 +271,7 @@ for jjj in range(3):
grid()
subplot(313, title="Acceleration along Z")
- scatter(pos[dist<dist_max], e_errz_s[dist<dist_max] , c=dist[dist<dist_max], label="Sorted", marker='o', s=1, linewidth=0, cmap='jet')
+ scatter(pos[ortho<ortho_max], e_errz_s[ortho<ortho_max] , c=ortho[ortho<ortho_max], label="Sorted", marker='o', s=1, linewidth=0, cmap='jet')
plot([-limit_exact, -limit_exact], [-plot_range, plot_range], 'k--')
plot([limit_exact, limit_exact], [-plot_range, plot_range], 'k--')
text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
@@ -191,21 +287,21 @@ for jjj in range(3):
figure(frameon=True)
subplot(311, title="Acceleration along X")
- scatter(dist, e_errx_s )
+ scatter(ortho, e_errx_s )
text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
#xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
#ylim(-0.01, 0.01)
grid()
subplot(312, title="Acceleration along Y")
- scatter(dist, e_erry_s )
+ scatter(ortho, e_erry_s )
text( 0., 0.005, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
#xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
#ylim(-0.01, 0.01)
grid()
subplot(313, title="Acceleration along Z")
- scatter(dist, e_errz_s , label="Sorted")
+ scatter(ortho, e_errz_s , label="Sorted")
text( 0., 0.005, "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")
#xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
@@ -256,12 +352,11 @@ for jjj in range(3):
# # 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.)]
-
-
+ 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")
diff --git a/examples/test_bh_sorted.c b/examples/test_bh_sorted.c
index de3928df60d0d8e19248955b174a74360357a066..8ec54aa138e9ff68b8177f8965087f2e92760bc4 100644
--- a/examples/test_bh_sorted.c
+++ b/examples/test_bh_sorted.c
@@ -38,11 +38,11 @@
/* Some local constants. */
#define cell_pool_grow 1000
-#define cell_maxparts 100
+#define cell_maxparts 1
#define task_limit 1e8
#define const_G 1 // 6.6738e-8
#define dist_min 0.5 /* Used for legacy walk only */
-#define dist_cutoff_ratio 1.5
+#define dist_cutoff_ratio 1.2
#define iact_pair_direct iact_pair_direct_sorted_multipole
#define ICHECK -1
#define NO_SANITY_CHECKS
@@ -60,6 +60,7 @@ struct part {
// };
float mass;
int id;
+ float d;
}; // __attribute__((aligned(64)));
struct part_local {
@@ -1336,10 +1337,12 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
for (i = 0; i < count_i; i++) {
int pid = ind_i[i].ind;
for (k = 0; k < 3; k++) ci->parts[pid].a[k] += parts_i[i].a[k];
+ ci->parts[pid].d = parts_i[i].d;
}
for (j = 0; j < count_j; j++) {
int pjd = ind_j[j].ind;
for (k = 0; k < 3; k++) cj->parts[pjd].a[k] += parts_j[j].a[k];
+ cj->parts[pjd].d = parts_j[j].d;
}
}
@@ -2083,13 +2086,13 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
file = fopen("particle_dump.dat", "w");
fprintf(file,
"# ID m x y z a_exact.x a_exact.y a_exact.z a_legacy.x "
- "a_legacy.y a_legacy.z a_new.x a_new.y a_new.z\n");
+ "a_legacy.y a_legacy.z a_new.x a_new.y a_new.z d\n");
for (k = 0; k < N; ++k)
- fprintf(file, "%d %e %e %e %e %e %e %e %e %e %e %e %e %e\n", parts[k].id,
+ fprintf(file, "%d %e %e %e %e %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_legacy[0], parts[k].a_legacy[1], parts[k].a_legacy[2],
- parts[k].a[0], parts[k].a[1], parts[k].a[2]);
+ parts[k].a[0], parts[k].a[1], parts[k].a[2], parts[k].d);
fclose(file);
/* Clean up. */
@@ -2168,7 +2171,7 @@ void test_direct_neighbour(int N_parts, int orientation) {
parts[k].x[2] += shift[2];
}
- parts[k].mass = ((double)rand()) / RAND_MAX;
+ parts[k].mass = 1.;//((double)rand()) / RAND_MAX;
parts[k].a[0] = 0.0;
parts[k].a[1] = 0.0;
parts[k].a[2] = 0.0;
@@ -2278,13 +2281,13 @@ void test_direct_neighbour(int N_parts, int orientation) {
message("Writing file '%s'", fileName);
FILE *file = fopen(fileName, "w");
fprintf(file,
- "# ID m r x y z a_u.x a_u.y a_u.z a_s.x a_s.y a_s.z ortho\n");
+ "# ID m r x y z a_u.x a_u.y a_u.z a_s.x a_s.y a_s.z ortho d\n");
for (k = 0; k < 2 * N_parts; k++) {
- fprintf(file, "%d %e %e %e %e %e %e %e %e %e %e %e %e\n", parts[k].id,
+ fprintf(file, "%d %e %e %e %e %e %e %e %e %e %e %e %e %e\n", parts[k].id,
parts[k].mass, position[k], 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],
- ortho_dist[k]);
+ ortho_dist[k], parts[k].d);
}
fclose(file);
@@ -2350,7 +2353,7 @@ void test_all_direct_neighbours(int N_parts, int N_test) {
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].mass = 1.;//((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;
@@ -2419,7 +2422,7 @@ void test_all_direct_neighbours(int N_parts, int N_test) {
qsched_init(&s, 1, qsched_flag_noreown);
cell_split(root, &s);
legacy_tree_walk(27 * N_parts, parts, root, ICHECK, &countMultipoles,
- &countPairs, &countCoMs);
+ &countPairs, &countCoMs);
// free(root);
//++countCoMs;