diff --git a/examples/plot_sorted.py b/examples/plot_sorted.py
index 9e07de4edc6d0afad13067742e1a929ed21c8b05..91de1340e04b155d57ad12242eec38b368ad1926 100644
--- a/examples/plot_sorted.py
+++ b/examples/plot_sorted.py
@@ -79,15 +79,14 @@ limit_exact = ( dist_cutoff_ratio - 1. )
print limit_exact
#names = ["side", "edge", "corner"]
-#for orientation in range( 26 ):
-for jjj in range(0,3):
- if jjj == 0:
- orientation = 0
- if jjj == 1:
- orientation = 1
- if jjj == 2:
- orientation = 4
-
+for orientation in range( 26 ):
+# for jjj in range(0,3):
+# if jjj == 0:
+# orientation = 0
+# if jjj == 1:
+# orientation = 1
+# if jjj == 2:
+# orientation = 4
# Read Quickshed accelerations
data=loadtxt( "interaction_dump_%d.dat"%orientation )
@@ -101,7 +100,7 @@ for jjj in range(0,3):
ortho = data[:,12]
dist = data[:,13]
-
+
accx_u=data[:,6]
accy_u=data[:,7]
accz_u=data[:,8]
@@ -128,10 +127,28 @@ for jjj in range(0,3):
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
-
+ this_axis = array(axis[orientation*3:orientation*3+3])
+ #if this_axis[1] != 0.:
+ # continue
+ cell_edge = sqrt(sum(abs(array(this_axis))))# - sum(abs(this_axis[this_axis<0]))
+ print this_axis
+ print 2*cell_edge
+ print "actual:", orientation, min(dist), max(dist)
+ print "---"
+ continue
+
+ if cell_edge > 1.5:
+ cubic_param1 = [ 0.04730586, -0.37960746, 0.99371996, 0.15444646]
+ cubic_param2 = [-0.04397042, 0.10367151, -0.06130978, 1.00754152]
+
+ elif cell_edge > 1.1:
+ cubic_param1 = [ 0.00536257, -0.06370355, 0.19954405, 0.81828777]
+ cubic_param2 = [-0.00677179, -0.01237365, 0.02703505, 0.99483005]
+
+ else:
+ cubic_param1 = [0.01082535, -0.0262587, 0.00342499, 1.0151183]
+ cubic_param2 = [-0.00352327, 0.02770999, -0.02312047, 1.00282468]
+
figure()
xx = dist[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
@@ -147,22 +164,23 @@ for jjj in range(0,3):
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
+ # 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')
+ text( 0, 1.01, "axis=( %d %d %d )"%(axis[orientation*3 + 0], axis[orientation*3 + 1], axis[orientation*3 + 2]) , ha='center', backgroundcolor='w', fontsize=14)
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, 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-.')
+ #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)]
@@ -178,30 +196,32 @@ for jjj in range(0,3):
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
+ # 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, 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-.')
+ #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)
+ close()
+
# 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])
+ accx_s[i] /= cubic(dist[i], cubic_param2[0], cubic_param2[1], cubic_param2[2], cubic_param2[3])
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])
+ accx_s[i] /= cubic(dist[i], cubic_param1[0], cubic_param1[1], cubic_param1[2], cubic_param1[3])
# Build error ------------------------------------------------
@@ -284,32 +304,32 @@ for jjj in range(0,3):
- figure(frameon=True)
+ # figure(frameon=True)
- subplot(311, title="Acceleration along X")
- 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(311, title="Acceleration along X")
+ # 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(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(312, title="Acceleration along Y")
+ # 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(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)))
- #ylim(-0.01, 0.01)
- grid()
+ # subplot(313, title="Acceleration along Z")
+ # 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)))
+ # #ylim(-0.01, 0.01)
+ # grid()
- savefig("error_distance_%d.png"%orientation)
- close()
+ # savefig("error_distance_%d.png"%orientation)
+ # close()
diff --git a/examples/test_bh_sorted.c b/examples/test_bh_sorted.c
index 4b4c4edfda21a2b2fe6f7b3e094865b89583e15a..1ee7b327bad5558b3dbf67dc95166b8c1aa27f15 100644
--- a/examples/test_bh_sorted.c
+++ b/examples/test_bh_sorted.c
@@ -281,6 +281,69 @@ const int axis_sid[27] = {
/* ( 1 , 1 , 1 ) */ 0
};
+const float axis_min_dist[27] = {
+ /* ( -1 , -1 , -1 ) */ 0.f,
+ /* ( -1 , -1 , 0 ) */ 0.f,
+ /* ( -1 , -1 , 1 ) */ -5.773502691896258e-01,
+ /* ( -1 , 0 , -1 ) */ 0.f,
+ /* ( -1 , 0 , 0 ) */ 0.f,
+ /* ( -1 , 0 , 1 ) */ -7.071067811865475e-01f,
+ /* ( -1 , 1 , -1 ) */ -5.773502691896258e-01f,
+ /* ( -1 , 1 , 0 ) */ -7.071067811865475e-01f,
+ /* ( -1 , 1 , 1 ) */ -1.1547005383792516f,
+ /* ( 0 , -1 , -1 ) */ 0.f,
+ /* ( 0 , -1 , 0 ) */ 0.f,
+ /* ( 0 , -1 , 1 ) */ -7.071067811865475e-01f,
+ /* ( 0 , 0 , -1 ) */ 0.f,
+ /* ( 0 , 0 , 0 ) */ 0, /* Should never happen */
+ /* ( 0 , 0 , 1 ) */ -1.73205080756887729352f,
+ /* ( 0 , 1 , -1 ) */ -1.41421356237309504880f,
+ /* ( 0 , 1 , 0 ) */ -2.30940107675850305803f,
+ /* ( 0 , 1 , 1 ) */ -1.41421356237309504880f,
+ /* ( 1 , -1 , -1 ) */ -1.f,
+ /* ( 1 , -1 , 0 ) */ -2.12132034355964257320f,
+ /* ( 1 , -1 , 1 ) */ -2.30940107675850305803f,
+ /* ( 1 , 0 , -1 ) */ -2.12132034355964257320f,
+ /* ( 1 , 0 , 0 ) */ -2.88675134594812882254f,
+ /* ( 1 , 0 , 1 ) */ -1.41421356237309504880f,
+ /* ( 1 , 1 , -1 ) */ -1.f,
+ /* ( 1 , 1 , 0 ) */ -2.12132034355964257320f,
+ /* ( 1 , 1 , 1 ) */ -1.f
+
+
+};
+
+const float axis_max_dist[27] = {
+ /* ( -1 , -1 , -1 ) */ 3.46410161513775458704f,
+ /* ( -1 , -1 , 0 ) */ 2.82842712474619009760f,
+ /* ( -1 , -1 , 1 ) */ 2.88675134594812882254f,
+ /* ( -1 , 0 , -1 ) */ 2.82842712474619009760f,
+ /* ( -1 , 0 , 0 ) */ 2.f,
+ /* ( -1 , 0 , 1 ) */ 2.12132034355964257320f,
+ /* ( -1 , 1 , -1 ) */ 2.88675134594812882254f,
+ /* ( -1 , 1 , 0 ) */ 2.12132034355964257320f,
+ /* ( -1 , 1 , 1 ) */ 2.30940107675850305803f,
+ /* ( 0 , -1 , -1 ) */ 2.82842712474619009760f,
+ /* ( 0 , -1 , 0 ) */ 2.f,
+ /* ( 0 , -1 , 1 ) */ 2.12132034355964257320f,
+ /* ( 0 , 0 , -1 ) */ 2.f,
+ /* ( 0 , 0 , 0 ) */ 0, /* Should never happen */
+ /* ( 0 , 0 , 1 ) */ 1.73205080756887729352f,
+ /* ( 0 , 1 , -1 ) */ 1.41421356237309504880f,
+ /* ( 0 , 1 , 0 ) */ 1.15470053837925152901f,
+ /* ( 0 , 1 , 1 ) */ 1.41421356237309504880f,
+ /* ( 1 , -1 , -1 ) */ 1.f,
+ /* ( 1 , -1 , 0 ) */ 0.70710678118654752440f,
+ /* ( 1 , -1 , 1 ) */ 1.15470053837925152901f,
+ /* ( 1 , 0 , -1 ) */ 0.70710678118654752440f,
+ /* ( 1 , 0 , 0 ) */ 0.57735026918962576450f,
+ /* ( 1 , 0 , 1 ) */ 1.41421356237309504880f,
+ /* ( 1 , 1 , -1 ) */ 1.f,
+ /* ( 1 , 1 , 0 ) */ 0.70710678118654752440f,
+ /* ( 1 , 1 , 1 ) */ 1.f
+
+ };
+
/* Some forward declarations. */
void comp_com(struct cell *c);
@@ -492,7 +555,7 @@ void cell_sort(struct cell *c, const float *axis, int aid) {
* @param ind_j Sorted indices of the cell @c cj.
*/
void get_axis(struct cell **ci, struct cell **cj, struct index **ind_i,
- struct index **ind_j) {
+ struct index **ind_j, float *min_dist, float* max_dist) {
float axis[3];
float dx[3];
@@ -504,6 +567,10 @@ void get_axis(struct cell **ci, struct cell **cj, struct index **ind_i,
aid = 3 * aid + ((dx[k] < 0) ? 0 : (dx[k] > 0) ? 2 : 1);
}
+ /* Get the minimal and maximal distance */
+ *min_dist = axis_min_dist[aid > 12 ? 26-aid: aid+14];
+ *max_dist = axis_max_dist[aid > 12 ? 26-aid: aid+14];
+
/* Flip the cells? */
if (axis_flip[aid]) {
struct cell *temp = *ci;
@@ -525,7 +592,6 @@ void get_axis(struct cell **ci, struct cell **cj, struct index **ind_i,
*ind_i = &(*ci)->indices[aid * ((*ci)->count + 1)];
*ind_j = &(*cj)->indices[aid * ((*cj)->count + 1)];
-
/* Make sure the sorts are ok. */
/* for (int k = 1; k < (*ci)->count; k++)
if ((*ind_i)[k].d < (*ind_i)[k-1].d)
@@ -1146,7 +1212,15 @@ static inline void iact_pair_direct_unsorted(struct cell *ci, struct cell *cj) {
-
+float correction_coefs[6*4] =
+ {
+ 0.04730586, -0.37960746, 0.99371996, 0.15444646,
+ -0.04397042, 0.10367151, -0.06130978, 1.00754152,
+ 0.00536257, -0.06370355, 0.19954405, 0.81828777,
+ -0.00677179, -0.01237365, 0.02703505, 0.99483005,
+ 0.01082535, -0.0262587, 0.00342499, 1.0151183,
+ -0.00352327, 0.02770999, -0.02312047, 1.00282468
+ };
static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
@@ -1158,6 +1232,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
float cjh = cj->h;
float xi[3];
float dx[3], ai[3], mi, mj, r2, w, ir;
+ float min_dist, max_dist;
#ifdef SANITY_CHECKS
@@ -1170,7 +1245,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
/* Get the sorted indices and stuff. */
struct index *ind_i, *ind_j;
struct multipole multi;
- get_axis(&ci, &cj, &ind_i, &ind_j);
+ get_axis(&ci, &cj, &ind_i, &ind_j, &min_dist, &max_dist);
multipole_reset(&multi);
cjh = cj->h;
@@ -2546,10 +2621,10 @@ int main(int argc, char *argv[]) {
N_parts);
/* Run the test */
- //for (k = 0; k < 26; ++k) test_direct_neighbour(N_parts, k);
- test_direct_neighbour(N_parts, 0);
- test_direct_neighbour(N_parts, 1);
- test_direct_neighbour(N_parts, 4);
+ for (k = 0; k < 26; ++k) test_direct_neighbour(N_parts, k);
+ //test_direct_neighbour(N_parts, 0);
+ //test_direct_neighbour(N_parts, 1);
+ //test_direct_neighbour(N_parts, 4);
k++;
/* Dump arguments */