diff --git a/examples/test_bh_new.c b/examples/test_bh_new.c
index 89c8caf3c9f49f5631fef378b7fc669bc97cb366..0fabd158cbaaf8cd5e53a03e2967ab4d57103ef2 100644
--- a/examples/test_bh_new.c
+++ b/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)