Skip to content
Snippets Groups Projects
Commit 7f599d35 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

The position along the axis linking the particles is now computed correctly.

parent 92ea88d1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
examples/plot_sorted.py
+ 80
− 60
View file @ 7f599d35
...@@ -79,15 +79,14 @@ limit_exact = ( dist_cutoff_ratio - 1. ) ...@@ -79,15 +79,14 @@ limit_exact = ( dist_cutoff_ratio - 1. )
print limit_exact print limit_exact
#names = ["side", "edge", "corner"] #names = ["side", "edge", "corner"]
#for orientation in range( 26 ): for orientation in range( 26 ):
for jjj in range(0,3): # for jjj in range(0,3):
if jjj == 0: # if jjj == 0:
orientation = 0 # orientation = 0
if jjj == 1: # if jjj == 1:
orientation = 1 # orientation = 1
if jjj == 2: # if jjj == 2:
orientation = 4 # orientation = 4
# Read Quickshed accelerations # Read Quickshed accelerations
data=loadtxt( "interaction_dump_%d.dat"%orientation ) data=loadtxt( "interaction_dump_%d.dat"%orientation )
...@@ -101,7 +100,7 @@ for jjj in range(0,3): ...@@ -101,7 +100,7 @@ for jjj in range(0,3):
ortho = data[:,12] ortho = data[:,12]
dist = data[:,13] dist = data[:,13]
accx_u=data[:,6] accx_u=data[:,6]
accy_u=data[:,7] accy_u=data[:,7]
accz_u=data[:,8] accz_u=data[:,8]
...@@ -128,10 +127,28 @@ for jjj in range(0,3): ...@@ -128,10 +127,28 @@ for jjj in range(0,3):
def quintic(x,a,b,c,d,e,f): 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 return a*x**5 + b*x**4 + c*x**3 + d*x**2 + e*x + f
this_axis = array(axis[orientation*3:orientation*3+3])
cell_edge = sqrt(sum(abs(array(axis[orientation*3:orientation*3+3])))) #if this_axis[1] != 0.:
print cell_edge # 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() figure()
xx = dist[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)] xx = dist[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
...@@ -147,22 +164,23 @@ for jjj in range(0,3): ...@@ -147,22 +164,23 @@ for jjj in range(0,3):
xx = xx[abs(cell_edge-xx) > limit_exact] xx = xx[abs(cell_edge-xx) > limit_exact]
delta = sim / 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) # 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 # #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) # 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 # 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) # 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 # #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) # 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 # #print quintic_param1
plot(xx, delta, 'bx') plot(xx, delta, 'bx')
plot([0,2*cell_edge], [1,1], 'r') 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) 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, 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, 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, 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)] xx = dist[(ortho < ortho_max_fit) & (ortho > ortho_min_fit)]
...@@ -178,30 +196,32 @@ for jjj in range(0,3): ...@@ -178,30 +196,32 @@ for jjj in range(0,3):
xx = xx[abs(cell_edge-xx) > limit_exact] xx = xx[abs(cell_edge-xx) > limit_exact]
delta = sim / 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) # 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 # #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) # 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 # 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) # 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 # #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) # 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 # #print quintic_param2
plot(xx, delta, 'bx') plot(xx, delta, 'bx')
xxx = linspace(0., cell_edge-limit_exact) 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, 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, 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, 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) savefig("fit_%d.png"%orientation)
close()
# Apply the correction --------------------------------------- # Apply the correction ---------------------------------------
for i in range(size(x)): for i in range(size(x)):
if cell_edge-dist[i] > limit_exact: 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: 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 ------------------------------------------------ # Build error ------------------------------------------------
...@@ -284,32 +304,32 @@ for jjj in range(0,3): ...@@ -284,32 +304,32 @@ for jjj in range(0,3):
figure(frameon=True) # figure(frameon=True)
subplot(311, title="Acceleration along X") # subplot(311, title="Acceleration along X")
scatter(ortho, 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) # 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))) # #xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
#ylim(-0.01, 0.01) # #ylim(-0.01, 0.01)
grid() # grid()
subplot(312, title="Acceleration along Y") # subplot(312, title="Acceleration along Y")
scatter(ortho, 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) # 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))) # #xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
#ylim(-0.01, 0.01) # #ylim(-0.01, 0.01)
grid() # grid()
subplot(313, title="Acceleration along Z") # subplot(313, title="Acceleration along Z")
scatter(ortho, 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) # 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") # #legend(loc="upper right")
#xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos))) # #xlim(-1.2*max(abs(pos)), 1.2*max(abs(pos)))
#ylim(-0.01, 0.01) # #ylim(-0.01, 0.01)
grid() # grid()
savefig("error_distance_%d.png"%orientation) # savefig("error_distance_%d.png"%orientation)
close() # close()
......
examples/test_bh_sorted.c
+ 83
− 8
View file @ 7f599d35
...@@ -281,6 +281,69 @@ const int axis_sid[27] = { ...@@ -281,6 +281,69 @@ const int axis_sid[27] = {
/* ( 1 , 1 , 1 ) */ 0 /* ( 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. */ /* Some forward declarations. */
void comp_com(struct cell *c); void comp_com(struct cell *c);
...@@ -492,7 +555,7 @@ void cell_sort(struct cell *c, const float *axis, int aid) { ...@@ -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. * @param ind_j Sorted indices of the cell @c cj.
*/ */
void get_axis(struct cell **ci, struct cell **cj, struct index **ind_i, 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 axis[3];
float dx[3]; float dx[3];
...@@ -504,6 +567,10 @@ void get_axis(struct cell **ci, struct cell **cj, struct index **ind_i, ...@@ -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); 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? */ /* Flip the cells? */
if (axis_flip[aid]) { if (axis_flip[aid]) {
struct cell *temp = *ci; struct cell *temp = *ci;
...@@ -525,7 +592,6 @@ void get_axis(struct cell **ci, struct cell **cj, struct index **ind_i, ...@@ -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_i = &(*ci)->indices[aid * ((*ci)->count + 1)];
*ind_j = &(*cj)->indices[aid * ((*cj)->count + 1)]; *ind_j = &(*cj)->indices[aid * ((*cj)->count + 1)];
/* Make sure the sorts are ok. */ /* Make sure the sorts are ok. */
/* for (int k = 1; k < (*ci)->count; k++) /* for (int k = 1; k < (*ci)->count; k++)
if ((*ind_i)[k].d < (*ind_i)[k-1].d) 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) { ...@@ -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, 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, ...@@ -1158,6 +1232,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
float cjh = cj->h; float cjh = cj->h;
float xi[3]; float xi[3];
float dx[3], ai[3], mi, mj, r2, w, ir; float dx[3], ai[3], mi, mj, r2, w, ir;
float min_dist, max_dist;
#ifdef SANITY_CHECKS #ifdef SANITY_CHECKS
...@@ -1170,7 +1245,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci, ...@@ -1170,7 +1245,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
/* Get the sorted indices and stuff. */ /* Get the sorted indices and stuff. */
struct index *ind_i, *ind_j; struct index *ind_i, *ind_j;
struct multipole multi; 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); multipole_reset(&multi);
cjh = cj->h; cjh = cj->h;
...@@ -2546,10 +2621,10 @@ int main(int argc, char *argv[]) { ...@@ -2546,10 +2621,10 @@ int main(int argc, char *argv[]) {
N_parts); N_parts);
/* Run the test */ /* Run the test */
//for (k = 0; k < 26; ++k) test_direct_neighbour(N_parts, k); for (k = 0; k < 26; ++k) test_direct_neighbour(N_parts, k);
test_direct_neighbour(N_parts, 0); //test_direct_neighbour(N_parts, 0);
test_direct_neighbour(N_parts, 1); //test_direct_neighbour(N_parts, 1);
test_direct_neighbour(N_parts, 4); //test_direct_neighbour(N_parts, 4);
k++; k++;
/* Dump arguments */ /* Dump arguments */
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment