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

Unbiassed correction is now applied on all 3 axis in all 26 configurations.

parent bd32afdd
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
examples/plot_sorted.py
+ 39
− 34
View file @ b598c747
......@@ -22,7 +22,7 @@ params = {'axes.labelsize': 16,
'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.markersize' : 4,
'lines.linewidth' : 2,
# 'axes.formatter.limits' : (-2, 1),
......@@ -79,14 +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 )
......@@ -157,30 +157,30 @@ for jjj in range(0,3):
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[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]
sim = sim[xx > limit_exact]
exact = exact[xx > limit_exact]
xx = xx[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
#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')
plot([min(dist), max(dist)], [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(limit_exact, max(dist))
#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, quintic(xxx, quintic_param1[0], quintic_param1[1], quintic_param1[2], quintic_param1[3], quintic_param1[4], quintic_param1[5]), 'k-.')
......@@ -189,30 +189,35 @@ for jjj in range(0,3):
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[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]
sim = sim[xx < -limit_exact]
exact = exact[xx < -limit_exact]
xx = xx[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
#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)
xxx = linspace(min(dist), -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, 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([-limit_exact, -limit_exact], [0.9, 1.1], 'k--')
plot([limit_exact, limit_exact], [0.9, 1.1], 'k--')
ylim(0.99, 1.01)
savefig("fit_%d.png"%orientation)
close()
......@@ -220,12 +225,12 @@ for jjj in range(0,3):
# Apply the correction ---------------------------------------
# for i in range(size(x)):
# if cell_edge-dist[i] > limit_exact:
# if dist[i] < -limit_exact:
# 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] > limit_exact:
# accx_s[i] /= cubic(dist[i], cubic_param1[0], cubic_param1[1], cubic_param1[2], cubic_param1[3])
# Build error ------------------------------------------------
inv_acc_tot = 1.0 / sqrt(accx_u**2 + accy_u**2 + accz_u**2)
......
examples/test_bh_sorted.c
+ 38
− 18
View file @ b598c747
......@@ -1166,12 +1166,12 @@ 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.00701166, -0.0113404, -0.0159732, 1.02330933,
0.00297062, 0.01928414, -0.01929454, 1.00190972
0.04730587, -0.13379902, 0.10447393, 0.982606,
-0.04397043, -0.12480554, -0.09791498, 0.98388736,
0.00536261, -0.04095222, 0.05153867, 0.98824618,
-0.00677177, -0.04110387, -0.04859363, 0.98916257,
0.00425108, 0.01519535, -0.02025942, 1.0035729,
-0.00711434, 0.01114372, 0.01865343, 1.00334988
};
......@@ -1216,16 +1216,16 @@ static inline void gravity_bias_correction( float d, float axis[3], float corr[3
case 1: /* side */
{
if ( d > limit_exact)
for( k = 0; k < 3; k++ )
corr[k] = 1./cubic( d, correction_coefs[4*4+0], correction_coefs[4*4+1], correction_coefs[4*4+2], correction_coefs[4*4+3]);
else if ( d < -limit_exact)
if ( d < -limit_exact )
for( k = 0; k < 3; k++ )
corr[k] = 1./cubic( d, correction_coefs[5*4+0], correction_coefs[5*4+1], correction_coefs[5*4+2], correction_coefs[5*4+3]);
corr[k] = 1./cubic(d, correction_coefs[5*4+0], correction_coefs[5*4+1], correction_coefs[5*4+2], correction_coefs[5*4+3]);
message("corr= %f %f %f", corr[0], corr[1], corr[2]);
else if (d > limit_exact)
for( k = 0; k < 3; k++ )
corr[k] = 1./cubic(d, correction_coefs[4*4+0], correction_coefs[4*4+1], correction_coefs[4*4+2], correction_coefs[4*4+3]);
}
break;
......@@ -1233,9 +1233,13 @@ static inline void gravity_bias_correction( float d, float axis[3], float corr[3
case 2: /* edge */
{
if (d > limit_exact || d < -limit_exact)
if ( d < -limit_exact )
for( k = 0; k < 3; k++ )
corr[k] = 1./cubic(d, correction_coefs[3*4+0], correction_coefs[3*4+1], correction_coefs[3*4+2], correction_coefs[3*4+3]);
else if (d > limit_exact)
for( k = 0; k < 3; k++ )
corr[k] = 0.;
corr[k] = 1./cubic(d, correction_coefs[2*4+0], correction_coefs[2*4+1], correction_coefs[2*4+2], correction_coefs[2*4+3]);
}
break;
......@@ -1243,9 +1247,14 @@ static inline void gravity_bias_correction( float d, float axis[3], float corr[3
case 3: /* corner */
{
if (d > limit_exact || d < -limit_exact)
if ( d < -limit_exact )
for( k = 0; k < 3; k++ )
corr[k] = 0.;
corr[k] = 1./cubic(d, correction_coefs[1*4+0], correction_coefs[1*4+1], correction_coefs[1*4+2], correction_coefs[1*4+3]);
else if (d > limit_exact)
for( k = 0; k < 3; k++ )
corr[k] = 1./cubic(d, correction_coefs[0*4+0], correction_coefs[0*4+1], correction_coefs[0*4+2], correction_coefs[0*4+3]);
}
break;
......@@ -1469,10 +1478,16 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
/* Get the correction */
gravity_bias_correction( d_norm, axis, corr );
//corr[0] = 1.;
//corr[1] = 1.;
//corr[2] = 1.;
/* And apply it */
parts_i[i].a[0] *= corr[0];
parts_i[i].a[1] *= corr[1];
parts_i[i].a[2] *= corr[2];
parts_i[i].d = d_norm;
}
......@@ -1486,10 +1501,15 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
/* Get the correction */
gravity_bias_correction( d_norm, axis, corr );
//corr[0] = 1.;
//corr[1] = 1.;
//corr[2] = 1.;
/* And apply it */
parts_j[j].a[0] *= corr[0];
parts_j[j].a[1] *= corr[1];
parts_j[j].a[2] *= corr[2];
parts_j[j].d = d_norm;
}
/* Copy the accelerations back to the original particles. */
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment