diff --git a/examples/plot_sorted.py b/examples/plot_sorted.py
index d66fec58514b706357492bb10d21d3222ad7a254..2caf184b597627575784e4bbfdf9a9101a9681a9 100644
--- a/examples/plot_sorted.py
+++ b/examples/plot_sorted.py
@@ -36,7 +36,7 @@ import os
from scipy import stats
from scipy import optimize
-dist_cutoff_ratio=1.2
+dist_cutoff_ratio=1.5
ortho_max_fit = 0.2
ortho_min_fit = 0.
ortho_max = 2.5
@@ -137,7 +137,7 @@ for orientation in range( 26 ):
print this_axis
print "actual:", orientation, min(dist), max(dist), min_dist[0], max_dist[0]
print "---"
- continue
+ #continue
if cell_edge > 1.5:
cubic_param1 = [ 0.04730586, -0.37960746, 0.99371996, 0.15444646]
@@ -219,11 +219,11 @@ for orientation in range( 26 ):
# Apply the correction ---------------------------------------
- for i in range(size(x)):
- if cell_edge-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:
- accx_s[i] /= cubic(dist[i], cubic_param1[0], cubic_param1[1], cubic_param1[2], cubic_param1[3])
+ # for i in range(size(x)):
+ # if cell_edge-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:
+ # accx_s[i] /= cubic(dist[i], cubic_param1[0], cubic_param1[1], cubic_param1[2], cubic_param1[3])
# Build error ------------------------------------------------
diff --git a/examples/test_bh_sorted.c b/examples/test_bh_sorted.c
index 895d079680bef6b752791f10bc74cc41d5667201..2345aa11f22bfea1a807f62f02777aeac67bc6de 100644
--- a/examples/test_bh_sorted.c
+++ b/examples/test_bh_sorted.c
@@ -42,7 +42,7 @@
#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.2
+#define dist_cutoff_ratio 1.5
#define iact_pair_direct iact_pair_direct_sorted_multipole
#define ICHECK -1
#define NO_SANITY_CHECKS
@@ -495,9 +495,8 @@ 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, float *min_dist, float* max_dist) {
+ struct index **ind_j, float *min_dist, float* max_dist, float axis[3]) {
- float axis[3];
float dx[3];
int aid = 0;
@@ -1176,6 +1175,78 @@ float correction_coefs[6*4] =
};
+/**
+ * @brief Compute the best-fitting bias correction factor for the acceleration of a particle
+ *
+ * @param d The normalized position of the particle along the axis
+ * @param axis The axis along which the two cells are interacting
+ * @param The correction to apply in each direction
+ */
+
+static inline void gravity_bias_correction( float d, float axis[3], float corr[3] ) {
+
+ int k, orientation;
+ int is_zero[3];
+ float limit_exact = ( dist_cutoff_ratio - 1. );
+
+ /* Start with no applied correction */
+ corr[0] = corr[1] = corr[2] = 1.f;
+
+ /* Find which case we are dealing with */
+ orientation = 0 ;
+ for (k = 0 ; k < 3; ++k) {
+ is_zero[k] = ( fabs(axis[k]) > 0 ? 1 : 0 );
+ orientation += is_zero[k];
+ }
+
+ //message("%f %f", d, limit_exact);
+
+ /* Apply the correction corresponding to this axis */
+ switch(orientation)
+ {
+
+ case 1: /* side */
+ {
+ if (d > limit_exact || d < -limit_exact)
+ for( k = 0; k < 3; k++ )
+ corr[k] = 0.;
+
+ }
+ break;
+
+ case 2: /* edge */
+ {
+
+ if (d > limit_exact || d < -limit_exact)
+ for( k = 0; k < 3; k++ )
+ corr[k] = 0.;
+
+ }
+ break;
+
+ case 3: /* corner */
+ {
+
+ if (d > limit_exact || d < -limit_exact)
+ for( k = 0; k < 3; k++ )
+ corr[k] = 0.;
+
+ }
+ break;
+
+ default:
+ error( "Invalid orientation axis=(%f %f %f) orientation= %d", axis[0], axis[1], axis[2], orientation);
+
+ }
+
+ /* Correction only applies in the direction along the axis. Rest is untouched */
+ for ( k = 0; k < 3; ++k )
+ corr[k] = (is_zero[k] ? corr[k] : 1 );
+
+
+}
+
+
static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
struct cell *cj) {
@@ -1184,7 +1255,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
struct part_local *parts_i, *parts_j;
float xi[3];
float dx[3], ai[3], mi, mj, r2, w, ir;
- float min_dist, max_dist, length_dist, mid_dist;
+ float min_dist, max_dist, length_dist, mid_dist, d_norm, axis[3], corr[3];
#ifdef SANITY_CHECKS
@@ -1196,7 +1267,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
/* Get the sorted indices and geometry. */
struct index *ind_i, *ind_j;
- get_axis(&ci, &cj, &ind_i, &ind_j, &min_dist, &max_dist);
+ get_axis(&ci, &cj, &ind_i, &ind_j, &min_dist, &max_dist, axis);
length_dist = max_dist - min_dist;
mid_dist = min_dist + 0.5f * length_dist;
@@ -1373,9 +1444,36 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
}
/* Apply the correction to cancel the bias */
-
for (i = 0; i < count_i; i++) {
- message("Length = %f Mid = %f dist = %f", length_dist, mid_dist, parts_i[i].d);
+
+ /* Get the normalized position along the axis */
+ d_norm = parts_i[i].d - mid_dist;
+ d_norm /= ci->h;
+
+ /* Get the correction */
+ gravity_bias_correction( d_norm, axis, corr );
+
+ /* And apply it */
+ parts_i[i].a[0] *= corr[0];
+ parts_i[i].a[1] *= corr[1];
+ parts_i[i].a[2] *= corr[2];
+
+ }
+
+ /* Do the same on the other particles */
+ for (j = 0; j < count_i; j++) {
+
+ /* Get the normalized position along the axis */
+ d_norm = parts_j[j].d - mid_dist;
+ d_norm /= cj->h;
+
+ /* Get the correction */
+ gravity_bias_correction( d_norm, axis, corr );
+
+ /* And apply it */
+ parts_j[j].a[0] *= corr[0];
+ parts_j[j].a[1] *= corr[1];
+ parts_j[j].a[2] *= corr[2];
}
/* Copy the accelerations back to the original particles. */