New version of the B-H code that uses quadrupoles for increased accuracy

examples/test_bh_new.c

@@ -37,16 +37,17 @@
 
 /* Some local constants. */
 #define cell_pool_grow 1000
-#define cell_maxparts 100
-#define task_limit 1e8
+#define cell_maxparts 128
+#define task_limit 1
 #define const_G 1    // 6.6738e-8
 #define dist_min 0.5 /* Used for legacy walk only */
 #define dist_cutoff_ratio 1.5
 
 #define ICHECK -1
-#define NO_SANITY_CHECKS
+#define SANITY_CHECKS
 #define NO_COM_AS_TASK
 #define NO_COUNTERS
+#define QUADRUPOLES
 
 /** Data structure for the particles. */
 struct part {
@@ -67,11 +68,21 @@ struct part_local {
 } __attribute__((aligned(32)));
 
 struct multipole {
+
+#ifdef QUADRUPOLES
+  double I[3][3];
+#endif
+
   double com[3];
   float mass;
 };
 
 struct multipole_local {
+
+#ifdef QUADRUPOLES
+  float I[3][3];
+#endif
+
   float com[3];
   float mass;
 };
@@ -103,7 +114,7 @@ struct cell {
   int res, com_tid;
   struct index *indices;
 
-} __attribute__((aligned(128)));
+};// __attribute__((aligned(256)));
 
 /** Task types. */
 enum task_type {
@@ -173,6 +184,15 @@ void comp_com(struct cell *c) {
   struct part *parts = c->parts;
   double com[3] = {0.0, 0.0, 0.0}, mass = 0.0;
 
+#ifdef QUADRUPOLES
+  double I[3][3] = {{0.0, 0.0, 0.0} , {0.0, 0.0, 0.0} , {0.0, 0.0, 0.0} };
+#endif
+
+#ifdef SANITY_CHECKS
+  if ( count == 0 )
+    error("CoM computed for an empty cell");
+#endif
+  
   if (c->split) {
 
     /* Loop over the projenitors and collect the multipole information. */
@@ -182,9 +202,29 @@ void comp_com(struct cell *c) {
       com[1] += cp->new.com[1] * cp_mass;
       com[2] += cp->new.com[2] * cp_mass;
       mass += cp_mass;
+
+#ifdef QUADRUPOLES
+      I[0][0] += cp->new.I[0][0];
+      I[0][1] += cp->new.I[0][1];
+      I[0][2] += cp->new.I[0][2];
+      I[1][0] += cp->new.I[1][0];
+      I[1][1] += cp->new.I[1][1];
+      I[1][2] += cp->new.I[1][2];
+      I[2][0] += cp->new.I[2][0];
+      I[2][1] += cp->new.I[2][1];
+      I[2][2] += cp->new.I[2][2];
+#endif 
+
     }
 
-    /* Otherwise, collect the multipole from the particles. */
+    /* Final operation for the COM */
+    float imass = 1.0f / mass;
+    com[0] *= imass;
+    com[1] *= imass;
+    com[2] *= imass;
+
+
+  /* Otherwise, collect the multipole from the particles. */
   } else {
 
     for (k = 0; k < count; k++) {
@@ -193,22 +233,47 @@ void comp_com(struct cell *c) {
       com[1] += parts[k].x[1] * p_mass;
       com[2] += parts[k].x[2] * p_mass;
       mass += p_mass;
+
+#ifdef QUADRUPOLES
+      I[0][0] += parts[k].x[0] * parts[k].x[0] * p_mass;
+      I[0][1] += parts[k].x[0] * parts[k].x[1] * p_mass;
+      I[0][2] += parts[k].x[0] * parts[k].x[2] * p_mass;
+      I[1][0] += parts[k].x[1] * parts[k].x[0] * p_mass;
+      I[1][1] += parts[k].x[1] * parts[k].x[1] * p_mass;
+      I[1][2] += parts[k].x[1] * parts[k].x[2] * p_mass;
+      I[2][0] += parts[k].x[2] * parts[k].x[0] * p_mass;
+      I[2][1] += parts[k].x[2] * parts[k].x[1] * p_mass;
+      I[2][2] += parts[k].x[2] * parts[k].x[2] * p_mass;
+#endif 
+
     }
-  }
 
-  /* Store the COM data, if any was collected. */
-  if (mass > 0.0) {
+    /* Final operation for the COM */
     float imass = 1.0f / mass;
-    c->new.com[0] = com[0] * imass;
-    c->new.com[1] = com[1] * imass;
-    c->new.com[2] = com[2] * imass;
-    c->new.mass = mass;
-  } else {
-    c->new.com[0] = 0.0;
-    c->new.com[1] = 0.0;
-    c->new.com[2] = 0.0;
-    c->new.mass = 0.0;
+    com[0] *= imass;
+    com[1] *= imass;
+    com[2] *= imass;
   }
+
+  /* Store the multipole data */
+  c->new.mass = mass;
+  c->new.com[0] = com[0];
+  c->new.com[1] = com[1];
+  c->new.com[2] = com[2];
+
+
+#ifdef QUADRUPOLES
+  c->new.I[0][0] = I[0][0];
+  c->new.I[0][1] = I[0][1];
+  c->new.I[0][2] = I[0][2];
+  c->new.I[1][0] = I[1][0];
+  c->new.I[1][1] = I[1][1];
+  c->new.I[1][2] = I[1][2];
+  c->new.I[2][0] = I[2][0];
+  c->new.I[2][1] = I[2][1];
+  c->new.I[2][2] = I[2][2];
+#endif
+
 }
 
 /**
@@ -427,7 +492,7 @@ static inline void make_interact_pc(struct part_local *parts, int part_count,
 				    struct multipole_local mults[8], int mult_count) {
 
   int i, j, k;
-  float dx[3] = {0.0, 0.0, 0.0}, r2, ir, w;
+  float dx[3] = {0.0, 0.0, 0.0}, r2, ir, mrinv3;
 
 #ifdef SANITY_CHECKS
 
@@ -455,8 +520,37 @@ static inline void make_interact_pc(struct part_local *parts, int part_count,
 
       /* Apply the gravitational acceleration. */
       ir = 1.0f / sqrtf(r2);
-      w = mults[j].mass * const_G * ir * ir * ir;
-      for (k = 0; k < 3; k++) parts[i].a[k] += w * dx[k];
+      mrinv3 = mults[j].mass * const_G * ir * ir * ir;
+
+#ifndef QUADRUPOLES
+
+	parts[i].a[0] += mrinv3 * dx[0];
+	parts[i].a[1] += mrinv3 * dx[1];
+	parts[i].a[2] += mrinv3 * dx[2];
+
+#else
+      
+      float mrinv5 = mrinv3 * ir * ir;
+      float mrinv7 = mrinv5 * ir * ir;
+
+      float D1 = -mrinv3;
+      float D2 =  3.f * mrinv5;
+      float D3 = -15.f * mrinv7;
+
+      float q = mults[j].I[0][0] + mults[j].I[1][1] + mults[j].I[2][2];
+      float qRx = mults[j].I[0][0]*dx[0] + mults[j].I[0][1]*dx[1] + mults[j].I[0][2]*dx[2];
+      float qRy = mults[j].I[1][0]*dx[0] + mults[j].I[1][1]*dx[1] + mults[j].I[1][2]*dx[2];
+      float qRz = mults[j].I[2][0]*dx[0] + mults[j].I[2][1]*dx[1] + mults[j].I[2][2]*dx[2];
+      float qRR = qRx * dx[0] + qRy * dx[1] + qRz * dx[2];
+      float C = D1 + 0.5f * D2 * q + 0.5f * D3 * qRR;
+
+      parts[i].a[0] -= C * dx[0] + D2 * qRx;
+      parts[i].a[1] -= C * dx[1] + D2 * qRy;
+      parts[i].a[2] -= C * dx[2] + D2 * qRz;
+
+      /* message(" %f %f %f", C * dx[0] + D2 * qRx, C * dx[1] + D2 * qRy, C * dx[2] + D2 * qRz); */
+
+#endif
 
     } /* loop over every multipoles. */
   } /* loop over every particle. */
@@ -517,6 +611,21 @@ static inline void iact_pair_pc(struct cell *ci, struct cell *cj) {
 	mult_local[mult_count].com[2] = cps->new.com[2] - cp->loc[2];
 	mult_local[mult_count].mass = cps->new.mass;
 
+#ifdef QUADRUPOLES
+
+	/* Final operation to get the quadrupoles */
+	double imass = 1. / cps->new.mass;
+	mult_local[mult_count].I[0][0] = cps->new.I[0][0]*imass - cps->new.com[0] * cps->new.com[0];
+	mult_local[mult_count].I[0][1] = cps->new.I[0][1]*imass - cps->new.com[0] * cps->new.com[1];
+	mult_local[mult_count].I[0][2] = cps->new.I[0][2]*imass - cps->new.com[0] * cps->new.com[2];
+	mult_local[mult_count].I[1][0] = cps->new.I[1][0]*imass - cps->new.com[1] * cps->new.com[0];
+	mult_local[mult_count].I[1][1] = cps->new.I[1][1]*imass - cps->new.com[1] * cps->new.com[1];
+	mult_local[mult_count].I[1][2] = cps->new.I[1][2]*imass - cps->new.com[1] * cps->new.com[2];
+	mult_local[mult_count].I[2][0] = cps->new.I[2][0]*imass - cps->new.com[2] * cps->new.com[0];
+	mult_local[mult_count].I[2][1] = cps->new.I[2][1]*imass - cps->new.com[2] * cps->new.com[1];
+	mult_local[mult_count].I[2][2] = cps->new.I[2][2]*imass - cps->new.com[2] * cps->new.com[2];
+#endif	
+
 	mult_count++;
       }
     }
@@ -1172,13 +1281,14 @@ void interact_exact(int N, struct part *parts, int monitor) {
  */
 void test_bh(int N, int nr_threads, int runs, char *fileName) {
 
-  int i, k;
+  int i,j, k;
   struct cell *root;
   struct part *parts;
   FILE *file;
   struct qsched s;
   ticks tic, toc_run, tot_setup = 0, tot_run = 0;
   int countMultipoles = 0, countPairs = 0, countCoMs = 0;
+  //  double I[3][3] = {{0., 0., 0.},{0., 0., 0.},{0., 0., 0.}};
 
   /* Runner function. */
   void runner(int type, void * data) {
@@ -1231,6 +1341,23 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
       parts[k].a_legacy[1] = 0.0;
       parts[k].a_legacy[2] = 0.0;
     }
+    j++;
+    /* for (i = 0; i < 4; i++) { */
+    /*   for (j = 0; j < 4; j++) { */
+    /* 	for (k = 0; k < 4; k++) { */
+    /* 	  int l = 16*i + 4*j + k; */
+    /* 	  parts[l].id = l; */
+    /* 	  parts[l].x[0] = 0.125 + 0.25 * i;  */
+    /* 	  parts[l].x[1] = 0.125 + 0.25 * j; */
+    /* 	  parts[l].x[2] = 0.125 + 0.25 * k; */
+    /* 	  parts[l].mass = 1.; */
+    /* 	  parts[l].a_legacy[0] = 0.0; */
+    /* 	  parts[l].a_legacy[1] = 0.0; */
+    /* 	  parts[l].a_legacy[2] = 0.0; */
+
+    /* 	} */
+    /*   } */
+    /* } */
 
     /* Otherwise, read them from a file. */
   } else {
@@ -1274,7 +1401,7 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
   }
   message("Average number of parts per leaf is %f.", ((double)N) / nr_leaves);
 
-#if ICHECK > 0
+#if 1 //ICHECK > 0
   printf("----------------------------------------------------------\n");
 
   /* Do a N^2 interactions calculation */
@@ -1355,10 +1482,43 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
     tot_run += toc_run - tic;
   }
 
-  message("[check] root mass= %f %f", root->legacy.mass, root->new.mass);
-  message("[check] root CoMx= %f %f", root->legacy.com[0], root->new.com[0]);
-  message("[check] root CoMy= %f %f", root->legacy.com[1], root->new.com[1]);
-  message("[check] root CoMz= %f %f", root->legacy.com[2], root->new.com[2]);
+  /* for (k = 0; k < N; ++k) { */
+  /*   I[0][0] += parts[k].mass * parts[k].x[0] * parts[k].x[0]; */
+  /*   I[0][1] += parts[k].mass * parts[k].x[0] * parts[k].x[1]; */
+  /*   I[0][2] += parts[k].mass * parts[k].x[0] * parts[k].x[2]; */
+  /*   I[1][0] += parts[k].mass * parts[k].x[1] * parts[k].x[0]; */
+  /*   I[1][1] += parts[k].mass * parts[k].x[1] * parts[k].x[1]; */
+  /*   I[1][2] += parts[k].mass * parts[k].x[1] * parts[k].x[2]; */
+  /*   I[2][0] += parts[k].mass * parts[k].x[2] * parts[k].x[0]; */
+  /*   I[2][1] += parts[k].mass * parts[k].x[2] * parts[k].x[1];       */
+  /*   I[2][2] += parts[k].mass * parts[k].x[2] * parts[k].x[2]; */
+  /* } */
+
+  /* for (k = 0; k < N; ++k) { */
+  /*   I[0][0] = I[0][0] / root->new.mass - root->new.com[0]*root->new.com[0]; */
+  /*   I[0][1] = I[0][1] / root->new.mass - root->new.com[0]*root->new.com[1]; */
+  /*   I[0][2] = I[0][2] / root->new.mass - root->new.com[0]*root->new.com[2]; */
+  /*   I[1][0] = I[1][0] / root->new.mass - root->new.com[1]*root->new.com[0]; */
+  /*   I[1][1] = I[1][1] / root->new.mass - root->new.com[1]*root->new.com[1]; */
+  /*   I[1][2] = I[1][2] / root->new.mass - root->new.com[1]*root->new.com[2]; */
+  /*   I[2][0] = I[2][0] / root->new.mass - root->new.com[2]*root->new.com[0]; */
+  /*   I[2][1] = I[2][1] / root->new.mass - root->new.com[2]*root->new.com[1]; */
+  /*   I[2][2] = I[2][2] / root->new.mass - root->new.com[2]*root->new.com[2]; */
+      
+  /* } */
+
+  /* message("[check] root mass= %f %f", root->legacy.mass, root->new.mass); */
+  /* message("[check] root CoMx= %f %f", root->legacy.com[0], root->new.com[0]); */
+  /* message("[check] root CoMy= %f %f", root->legacy.com[1], root->new.com[1]); */
+  /* message("[check] root CoMz= %f %f", root->legacy.com[2], root->new.com[2]); */
+  /* message("   | %f %f %f |", root->new.I[0][0], root->new.I[0][1], root->new.I[0][2]); */
+  /* message("In=| %f %f %f |", root->new.I[1][0], root->new.I[1][1], root->new.I[1][2]); */
+  /* message("   | %f %f %f |", root->new.I[2][0], root->new.I[2][1], root->new.I[2][2]); */
+  /* message("------------------------------------------------") */
+  /* message("   | %f %f %f |", I[0][0], I[0][1], I[0][2]); */
+  /* message("I =| %f %f %f |", I[1][0], I[1][1], I[1][2]); */
+  /* message("   | %f %f %f |", I[2][0], I[2][1], I[2][2]); */
+
 #if ICHECK >= 0
   for (i = 0; i < N; ++i)
     if (root->parts[i].id == ICHECK)