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

Initial implementation of the FMM interaction and downward pass. Dowardpass...

Initial implementation of the FMM interaction and downward pass. Dowardpass needs to contain the acceleration shifts.
parent 091e277b
No related branches found
No related tags found
No related merge requests found
...@@ -58,7 +58,8 @@ struct part { ...@@ -58,7 +58,8 @@ struct part {
// }; // };
float mass; float mass;
int id; int id;
}; // __attribute__((aligned(32))); };// __attribute__((aligned(32)));
struct part_local { struct part_local {
float x[3]; float x[3];
...@@ -94,10 +95,12 @@ struct cell { ...@@ -94,10 +95,12 @@ struct cell {
struct multipole new; struct multipole new;
}; };
int res, com_tid; double a[3];
int res, com_tid, down_tid;
struct index *indices; struct index *indices;
} __attribute__((aligned(128))); } __attribute__((aligned(32)));
/** Task types. */ /** Task types. */
enum task_type { enum task_type {
...@@ -105,6 +108,7 @@ enum task_type { ...@@ -105,6 +108,7 @@ enum task_type {
task_type_pair, task_type_pair,
task_type_pair_pc, task_type_pair_pc,
task_type_com, task_type_com,
task_type_down,
task_type_count task_type_count
}; };
...@@ -203,8 +207,55 @@ void comp_com(struct cell *c) { ...@@ -203,8 +207,55 @@ void comp_com(struct cell *c) {
c->new.com[2] = 0.0; c->new.com[2] = 0.0;
c->new.mass = 0.0; c->new.mass = 0.0;
} }
/* Reset COM acceleration */
c->a[0] = 0.;
c->a[1] = 0.;
c->a[2] = 0.;
} }
/**
* @brief Add a cell's acceleration to all its siblings
*
* @param c The #cell.
*/
void comp_down(struct cell *c) {
int k, count = c->count;
struct cell *cp;
struct part *parts = c->parts;
// message("h=%f a= %f %f %f", c->h, c->a[0], c->a[1], c->a[2]);
if (c->split) {
//message("first= %p, sibling= %p", c->firstchild, c->sibling);
/* Loop over the siblings and propagate accelerations. */
for (cp = c->firstchild; cp != c->sibling; cp = cp->sibling) {
cp->a[0] += c->a[0];
cp->a[1] += c->a[1];
cp->a[2] += c->a[2];
/* Recurse */
//comp_down(cp);
}
} else {
/* Otherwise, propagate the accelerations to the siblings. */
for (k = 0; k < count; k++) {
parts[k].a[0] += c->a[0];
parts[k].a[1] += c->a[1];
parts[k].a[2] += c->a[2];
}
}
}
/** /**
* @brief Sort the parts into eight bins along the given pivots and * @brief Sort the parts into eight bins along the given pivots and
* fill the multipoles. Also adds the hierarchical resources * fill the multipoles. Also adds the hierarchical resources
...@@ -234,13 +285,19 @@ void cell_split(struct cell *c, struct qsched *s) { ...@@ -234,13 +285,19 @@ void cell_split(struct cell *c, struct qsched *s) {
if (c->res == qsched_res_none) if (c->res == qsched_res_none)
c->res = qsched_addres(s, qsched_owner_none, qsched_res_none); c->res = qsched_addres(s, qsched_owner_none, qsched_res_none);
/* Attach a center-of-mass task to the cell. */ /* Attach a center-of-mass and downward pass task to the cell. */
if (count > 0) {
#ifdef COM_AS_TASK #ifdef COM_AS_TASK
if (count > 0)
c->com_tid = qsched_addtask(s, task_type_com, task_flag_none, &c, c->com_tid = qsched_addtask(s, task_type_com, task_flag_none, &c,
sizeof(struct cell *), 1); sizeof(struct cell *), 1);
#endif #endif
c->down_tid = qsched_addtask(s, task_type_down, task_flag_none, &c,
sizeof(struct cell *), 1);
qsched_addlock(s, c->down_tid, c->res);
}
/* Does this cell need to be split? */ /* Does this cell need to be split? */
if (count > cell_maxparts) { if (count > cell_maxparts) {
...@@ -355,13 +412,21 @@ void cell_split(struct cell *c, struct qsched *s) { ...@@ -355,13 +412,21 @@ void cell_split(struct cell *c, struct qsched *s) {
for (k = 0; k < 8; k++) for (k = 0; k < 8; k++)
if (progenitors[k]->count > 0) cell_split(progenitors[k], s); if (progenitors[k]->count > 0) cell_split(progenitors[k], s);
/* Link the COM tasks. */ /* Link the COM tasks. */
#ifdef COM_AS_TASK #ifdef COM_AS_TASK
for (k = 0; k < 8; k++) for (k = 0; k < 8; k++)
if (progenitors[k]->count > 0) if (progenitors[k]->count > 0)
qsched_addunlock(s, progenitors[k]->com_tid, c->com_tid); qsched_addunlock(s, progenitors[k]->com_tid, c->com_tid);
#endif #endif
/* Link the downward tasks. */
#ifdef COM_AS_TASK
for (k = 0; k < 8; k++)
if (progenitors[k]->count > 0)
qsched_addunlock(s, c->down_tid, progenitors[k]->down_tid);
#endif
} /* does the cell need to be split? */ } /* does the cell need to be split? */
...@@ -408,58 +473,58 @@ static inline int are_neighbours(struct cell *ci, struct cell *cj) { ...@@ -408,58 +473,58 @@ static inline int are_neighbours(struct cell *ci, struct cell *cj) {
/** /* /\** */
* @brief Interacts all count particles in parts /* * @brief Interacts all count particles in parts */
* with the monopole in cj /* * with the monopole in cj */
* /* * */
* @param parts An array of particles /* * @param parts An array of particles */
* @param count The number of particles in the array /* * @param count The number of particles in the array */
* @param cj The cell with which the particles will interact /* * @param cj The cell with which the particles will interact */
*/ /* *\/ */
static inline void make_interact_pc(struct part *parts, int count, struct cell *cj) { /* static inline void make_interact_pc(struct part *parts, int count, struct cell *cj) { */
int j, k; /* int j, k; */
float com[3] = {0.0, 0.0, 0.0}; /* float com[3] = {0.0, 0.0, 0.0}; */
float mcom, dx[3] = {0.0, 0.0, 0.0}, r2, ir, w; /* float mcom, dx[3] = {0.0, 0.0, 0.0}, r2, ir, w; */
#ifdef SANITY_CHECKS /* #ifdef SANITY_CHECKS */
/* Sanity checks */ /* /\* Sanity checks *\/ */
if (count == 0) error("Empty cell!"); /* if (count == 0) error("Empty cell!"); */
/* Sanity check. */ /* /\* Sanity check. *\/ */
if (cj->new.mass == 0.0) { /* if (cj->new.mass == 0.0) { */
message("%e %e %e %d %p %d %p", cj->new.com[0], cj->new.com[1], /* message("%e %e %e %d %p %d %p", cj->new.com[0], cj->new.com[1], */
cj->new.com[2], cj->count, cj, cj->split, cj->sibling); /* cj->new.com[2], cj->count, cj, cj->split, cj->sibling); */
for (j = 0; j < cj->count; ++j) /* for (j = 0; j < cj->count; ++j) */
message("part %d mass=%e id=%d", j, cj->parts[j].mass, cj->parts[j].id); /* message("part %d mass=%e id=%d", j, cj->parts[j].mass, cj->parts[j].id); */
error("com does not seem to have been set."); /* error("com does not seem to have been set."); */
} /* } */
#endif /* #endif */
/* Init the com's data. */ /* /\* Init the com's data. *\/ */
for (k = 0; k < 3; k++) com[k] = cj->new.com[k]; /* for (k = 0; k < 3; k++) com[k] = cj->new.com[k]; */
mcom = cj->new.mass; /* mcom = cj->new.mass; */
/* Loop over every particle in leaf. */ /* /\* Loop over every particle in leaf. *\/ */
for (j = 0; j < count; j++) { /* for (j = 0; j < count; j++) { */
/* Compute the pairwise distance. */ /* /\* Compute the pairwise distance. *\/ */
for (r2 = 0.0, k = 0; k < 3; k++) { /* for (r2 = 0.0, k = 0; k < 3; k++) { */
dx[k] = com[k] - parts[j].x[k]; /* dx[k] = com[k] - parts[j].x[k]; */
r2 += dx[k] * dx[k]; /* r2 += dx[k] * dx[k]; */
} /* } */
/* Apply the gravitational acceleration. */ /* /\* Apply the gravitational acceleration. *\/ */
ir = 1.0f / sqrtf(r2); /* ir = 1.0f / sqrtf(r2); */
w = mcom * const_G * ir * ir * ir; /* w = mcom * const_G * ir * ir * ir; */
for (k = 0; k < 3; k++) parts[j].a[k] += w * dx[k]; /* for (k = 0; k < 3; k++) parts[j].a[k] += w * dx[k]; */
} /* loop over every particle. */ /* } /\* loop over every particle. *\/ */
} /* } */
...@@ -474,6 +539,8 @@ static inline void make_interact_pc(struct part *parts, int count, struct cell * ...@@ -474,6 +539,8 @@ static inline void make_interact_pc(struct part *parts, int count, struct cell *
static inline void iact_pair_pc(struct cell *ci, struct cell *cj) { static inline void iact_pair_pc(struct cell *ci, struct cell *cj) {
struct cell *cp, *cps; struct cell *cp, *cps;
int k;
float r2, w, ir, dx[3];
#ifdef SANITY_CHECKS #ifdef SANITY_CHECKS
if (ci->h != cj->h) if (ci->h != cj->h)
...@@ -489,8 +556,21 @@ static inline void iact_pair_pc(struct cell *ci, struct cell *cj) { ...@@ -489,8 +556,21 @@ static inline void iact_pair_pc(struct cell *ci, struct cell *cj) {
if ( ! are_neighbours(cp, cps) ) { if ( ! are_neighbours(cp, cps) ) {
make_interact_pc(cp->parts, cp->count, cps);
make_interact_pc(cps->parts, cps->count, cp );
/* Compute the pairwise distance. */
for (r2 = 0.0, k = 0; k < 3; k++) {
dx[k] = cps->new.com[k] - cp->new.com[k];
r2 += dx[k] * dx[k];
}
ir = 1.0f / sqrtf(r2);
/* Apply the gravitational acceleration. */
w = cps->new.mass * const_G * ir * ir * ir;
for (k = 0; k < 3; k++) cp->a[k] += w * dx[k];
w = cp->new.mass * const_G * ir * ir * ir;
for (k = 0; k < 3; k++) cps->a[k] -= w * dx[k];
} }
} }
...@@ -1155,6 +1235,9 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) { ...@@ -1155,6 +1235,9 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
case task_type_com: case task_type_com:
comp_com(d[0]); comp_com(d[0]);
break; break;
case task_type_down:
comp_down(d[0]);
break;
default: default:
error("Unknown task type."); error("Unknown task type.");
} }
...@@ -1227,7 +1310,7 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) { ...@@ -1227,7 +1310,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); message("Average number of parts per leaf is %f.", ((double)N) / nr_leaves);
#if 1//ICHECK > 0 #if 1 //ICHECK > 0
printf("----------------------------------------------------------\n"); printf("----------------------------------------------------------\n");
/* Do a N^2 interactions calculation */ /* Do a N^2 interactions calculation */
...@@ -1283,12 +1366,12 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) { ...@@ -1283,12 +1366,12 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
// fclose(fileTime); // fclose(fileTime);
printf("task counts: [ %10s %10s %10s %10s %10s ]\n", "self", "pair", "m-poles", printf("task counts: [ %10s %10s %10s %10s %10s %10s ]\n", "self", "pair", "m-poles",
"direct", "CoMs"); "direct", "CoMs", "Downpass");
printf("task counts: [ %10i %10i %10i %10i %10i ] (legacy).\n", 0, 0, printf("task counts: [ %10i %10i %10i %10i %10i %10i ] (legacy).\n", 0, 0,
countMultipoles, countPairs, countCoMs); countMultipoles, countPairs, countCoMs, 0);
printf("task counts: [ %10i %10i %10i %10i %10i ] (legacy).\n", counts[0], printf("task counts: [ %10i %10i %10i %10i %10i %10i ] (new).\n", counts[0],
counts[1], counts[2], 0, counts[3]); counts[1], counts[2], 0, counts[3], counts[4]);
/* Loop over the number of runs. */ /* Loop over the number of runs. */
for (k = 0; k < runs; k++) { for (k = 0; k < runs; k++) {
...@@ -1302,6 +1385,7 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) { ...@@ -1302,6 +1385,7 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
/* Execute the tasks. */ /* Execute the tasks. */
tic = getticks(); tic = getticks();
qsched_run(&s, nr_threads, runner); qsched_run(&s, nr_threads, runner);
comp_down(root);
toc_run = getticks(); toc_run = getticks();
message("%ith run took %lli (= %e) ticks...", k, toc_run - tic, message("%ith run took %lli (= %e) ticks...", k, toc_run - tic,
(float)(toc_run - tic)); (float)(toc_run - tic));
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment