Implemented 0-th order correction that just kills all inexact forces.

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 ------------------------------------------------

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. */