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

Back to a B-H implementation with loads of separated tasks, ready for FMM.

parent f8d4ff50
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 158 additions and 145 deletions
examples/test_fmm.c
+ 158
− 145
View file @ b6c23ed1
......@@ -38,7 +38,7 @@
/* Some local constants. */
#define cell_pool_grow 1000
#define cell_maxparts 100
#define task_limit 1e8
#define task_limit 1e10
#define const_G 1 // 6.6738e-8
#define dist_min 0.5 /* Used for legacy walk only */
#define dist_cutoff_ratio 1.5
......@@ -362,15 +362,17 @@ void cell_split(struct cell *c, struct qsched *s) {
qsched_addunlock(s, progenitors[k]->com_tid, c->com_tid);
#endif
/* Otherwise, we're at a leaf, so create the cell's particle-cell task. */
} else {
struct cell *data[2] = {root, c};
int tid = qsched_addtask(s, task_type_self_pc, task_flag_none, data,
2 * sizeof(struct cell *), 1);
qsched_addlock(s, tid, c->res);
#ifdef COM_AS_TASK
qsched_addunlock(s, root->com_tid, tid);
#endif
/* /\* Otherwise, we're at a leaf, so create the cell's particle-cell task. *\/ */
/* } else { */
/* struct cell *data[2] = {root, c}; */
/* int tid = qsched_addtask(s, task_type_self_pc, task_flag_none, data, */
/* 2 * sizeof(struct cell *), 1); */
/* qsched_addlock(s, tid, c->res); */
/* #ifdef COM_AS_TASK */
/* qsched_addunlock(s, root->com_tid, tid); */
/* #endif */
} /* does the cell need to be split? */
/* Compute the cell's center of mass. */
......@@ -390,7 +392,7 @@ void cell_split(struct cell *c, struct qsched *s) {
/**
* @brief Interacts all particles in ci with the monopole in cj
*/
static inline void make_interact_pc(struct part_local *parts, int count,
static inline void make_interact_pc(struct part *parts, int count,
double *loc, struct cell *cj) {
int j, k;
......@@ -416,7 +418,7 @@ static inline void make_interact_pc(struct part_local *parts, int count,
#endif
/* Init the com's data. */
for (k = 0; k < 3; k++) com[k] = cj->new.com[k] - loc[k];
for (k = 0; k < 3; k++) com[k] = cj->new.com[k];
mcom = cj->new.mass;
/* Loop over every particle in leaf. */
......@@ -488,148 +490,154 @@ static inline int are_neighbours(struct cell *ci, struct cell *cj) {
return 1;
}
/**
* @brief Compute the interactions between all particles in a cell leaf
* and the center of mass of all the cells in a part of the tree
* described by ci and cj
*
* @param ci The #cell containing the particle
* @param cj The #cell containing the center of mass.
*/
static inline void iact_pair_pc(struct cell *ci, struct cell *cj,
struct cell *leaf, struct part_local *parts,
int count, double *loc) {
/* /\** */
/* * @brief Compute the interactions between all particles in a cell leaf */
/* * and the center of mass of all the cells in a part of the tree */
/* * described by ci and cj */
/* * */
/* * @param ci The #cell containing the particle */
/* * @param cj The #cell containing the center of mass. */
/* *\/ */
/* static inline void iact_pair_pc(struct cell *ci, struct cell *cj, */
/* struct cell *leaf, struct part_local *parts, */
/* int count, double *loc) { */
struct cell *cp, *cps;
/* struct cell *cp, *cps; */
#ifdef SANITY_CHECKS
/* #ifdef SANITY_CHECKS */
/* Early abort? */
if (ci->count == 0 || cj->count == 0) error("Empty cell !");
/* /\* Early abort? *\/ */
/* if (ci->count == 0 || cj->count == 0) error("Empty cell !"); */
/* Sanity check */
if (ci == cj)
error("The impossible has happened: pair interaction between a cell and "
"itself.");
/* /\* Sanity check *\/ */
/* if (ci == cj) */
/* error("The impossible has happened: pair interaction between a cell and " */
/* "itself."); */
/* Sanity check */
if (!is_inside(leaf, ci))
error("The impossible has happened: The leaf is not within ci");
/* /\* Sanity check *\/ */
/* if (!is_inside(leaf, ci)) */
/* error("The impossible has happened: The leaf is not within ci"); */
/* Are the cells direct neighbours? */
if (!are_neighbours(ci, cj)) error("Cells are not neighours");
/* /\* Are the cells direct neighbours? *\/ */
/* if (!are_neighbours(ci, cj)) error("Cells are not neighours"); */
/* Are both cells split ? */
if (!ci->split || !cj->split) error("One of the cells is not split !");
#endif
/* /\* Are both cells split ? *\/ */
/* if (!ci->split || !cj->split) error("One of the cells is not split !"); */
/* #endif */
/* Let's find in which subcell of ci the leaf is */
for (cp = ci->firstchild; cp != ci->sibling; cp = cp->sibling) {
/* /\* Let's find in which subcell of ci the leaf is *\/ */
/* for (cp = ci->firstchild; cp != ci->sibling; cp = cp->sibling) { */
if (is_inside(leaf, cp)) break;
}
/* if (is_inside(leaf, cp)) break; */
/* } */
if (are_neighbours_different_size(cp, cj)) {
/* if (are_neighbours_different_size(cp, cj)) { */
/* Now interact this subcell with all subcells of cj */
for (cps = cj->firstchild; cps != cj->sibling; cps = cps->sibling) {
/* /\* Now interact this subcell with all subcells of cj *\/ */
/* for (cps = cj->firstchild; cps != cj->sibling; cps = cps->sibling) { */
/* Check whether we have to recurse or can directly jump to the multipole
* calculation */
if (are_neighbours(cp, cps)) {
/* /\* Check whether we have to recurse or can directly jump to the multipole */
/* * calculation *\/ */
/* if (are_neighbours(cp, cps)) { */
/* We only recurse if the children are split */
if (cp->split && cps->split) {
iact_pair_pc(cp, cps, leaf, parts, count, loc);
}
/* /\* We only recurse if the children are split *\/ */
/* if (cp->split && cps->split) { */
/* iact_pair_pc(cp, cps, leaf, parts, count, loc); */
/* } */
} else {
make_interact_pc(parts, count, loc, cps);
}
}
} else {
/* } else { */
/* make_interact_pc(parts, count, loc, cps); */
/* } */
/* } */
/* } else { */
/* If cp is not a neoghbour of cj, we can directly interact with the
* multipoles */
for (cps = cj->firstchild; cps != cj->sibling; cps = cps->sibling) {
/* /\* If cp is not a neoghbour of cj, we can directly interact with the */
/* * multipoles *\/ */
/* for (cps = cj->firstchild; cps != cj->sibling; cps = cps->sibling) { */
make_interact_pc(parts, count, loc, cps);
}
}
}
/* make_interact_pc(parts, count, loc, cps); */
/* } */
/* } */
/* } */
/**
* @brief Compute the interactions between all particles in a leaf and
* and all the monopoles in the cell c
*
* @param c The #cell containing the monopoles
* @param leaf The #cell containing the particles
*/
static inline void iact_self_pc(struct cell *c, struct cell *leaf,
struct part_local *parts) {
struct cell *cp, *cps;
int collect_part_data = 0;
#ifdef SANITY_CHECKS
/* Early abort? */
if (c->count == 0) error("Empty cell !");
/* /\** */
/* * @brief Compute the interactions between all particles in a leaf and */
/* * and all the monopoles in the cell c */
/* * */
/* * @param c The #cell containing the monopoles */
/* * @param leaf The #cell containing the particles */
/* *\/ */
/* static inline void iact_self_pc(struct cell *c, struct cell *leaf, */
/* struct part_local *parts) { */
if (!c->split) error("Cell is not split !");
/* struct cell *cp, *cps; */
/* int collect_part_data = 0; */
#endif
/* Get local copies of the particle data. */
if (parts == NULL) {
int count = leaf->count;
if ((parts =
(struct part_local *)malloc(sizeof(struct part_local) * count)) ==
NULL)
error("Failed to allocate local parts.");
for (int k = 0; k < count; k++) {
for (int j = 0; j < 3; j++) {
parts[k].x[j] = leaf->parts[k].x[j] - leaf->loc[j];
parts[k].a[j] = 0.0f;
}
parts[k].mass = leaf->parts[k].mass;
}
collect_part_data = 1;
}
/* #ifdef SANITY_CHECKS */
/* Find in which subcell of c the leaf is */
for (cp = c->firstchild; cp != c->sibling; cp = cp->sibling) {
/* /\* Early abort? *\/ */
/* if (c->count == 0) error("Empty cell !"); */
/* Only recurse if the leaf is in this part of the tree */
if (is_inside(leaf, cp)) break;
}
/* if (!c->split) error("Cell is not split !"); */
if (cp->split) {
/* Recurse if the cell can be split */
iact_self_pc(cp, leaf, parts);
/* Now, interact with every other subcell */
for (cps = c->firstchild; cps != c->sibling; cps = cps->sibling) {
/* #endif */
/* Since cp and cps will be direct neighbours it is only worth recursing
*/
/* if the cells can both be split */
if (cp != cps && cps->split)
iact_pair_pc(cp, cps, leaf, parts, leaf->count, leaf->loc);
}
}
/* /\* Get local copies of the particle data. *\/ */
/* if (parts == NULL) { */
/* int count = leaf->count; */
/* if ((parts = */
/* (struct part_local *)malloc(sizeof(struct part_local) * count)) == */
/* NULL) */
/* error("Failed to allocate local parts."); */
/* for (int k = 0; k < count; k++) { */
/* for (int j = 0; j < 3; j++) { */
/* parts[k].x[j] = leaf->parts[k].x[j] - leaf->loc[j]; */
/* parts[k].a[j] = 0.0f; */
/* } */
/* parts[k].mass = leaf->parts[k].mass; */
/* } */
/* collect_part_data = 1; */
/* } */
/* /\* Find in which subcell of c the leaf is *\/ */
/* for (cp = c->firstchild; cp != c->sibling; cp = cp->sibling) { */
/* /\* Only recurse if the leaf is in this part of the tree *\/ */
/* if (is_inside(leaf, cp)) break; */
/* } */
/* if (cp->split) { */
/* /\* Recurse if the cell can be split *\/ */
/* iact_self_pc(cp, leaf, parts); */
/* /\* Now, interact with every other subcell *\/ */
/* for (cps = c->firstchild; cps != c->sibling; cps = cps->sibling) { */
/* /\* Since cp and cps will be direct neighbours it is only worth recursing */
/* *\/ */
/* /\* if the cells can both be split *\/ */
/* if (cp != cps && cps->split) */
/* iact_pair_pc(cp, cps, leaf, parts, leaf->count, leaf->loc); */
/* } */
/* } */
/* /\* Clean up local parts? *\/ */
/* if (collect_part_data) { */
/* for (int k = 0; k < leaf->count; k++) { */
/* for (int j = 0; j < 3; j++) leaf->parts[k].a[j] += parts[k].a[j]; */
/* } */
/* free(parts); */
/* } */
/* } */
static inline void iact_pair_pc(struct cell *ci, struct cell *cj) {
make_interact_pc(ci->parts, ci->count, ci->loc, cj);
/* Clean up local parts? */
if (collect_part_data) {
for (int k = 0; k < leaf->count; k++) {
for (int j = 0; j < 3; j++) leaf->parts[k].a[j] += parts[k].a[j];
}
free(parts);
}
}
static inline void iact_pair_direct(struct cell *ci, struct cell *cj) {
int i, j, k;
......@@ -745,27 +753,32 @@ void iact_pair(struct cell *ci, struct cell *cj) {
error("The impossible has happened: pair interaction between a cell and "
"itself.");
#endif
/* Are the cells direct neighbours? */
if (are_neighbours(ci, cj)) {
if (!are_neighbours(ci, cj))
error("Non-neighbouring cells !");
/* Are both cells split ? */
if (ci->split && cj->split) {
#endif
/* Let's split both cells and build all possible pairs */
for (cp = ci->firstchild; cp != ci->sibling; cp = cp->sibling) {
for (cps = cj->firstchild; cps != cj->sibling; cps = cps->sibling) {
/* If the cells are neighbours, recurse. */
if (are_neighbours(cp, cps)) {
iact_pair(cp, cps);
}
}
/* Are both cells split ? */
if (ci->split && cj->split) {
/* Let's split both cells and build all possible pairs */
for (cp = ci->firstchild; cp != ci->sibling; cp = cp->sibling) {
for (cps = cj->firstchild; cps != cj->sibling; cps = cps->sibling) {
/* If the cells are neighbours, recurse. */
if (are_neighbours(cp, cps)) {
iact_pair(cp, cps);
}
else { /* Otherwise do cell-mulitpole interactions */
iact_pair_pc(cp, cps);
iact_pair_pc(cps, cp);
}
}
} else {/* Otherwise, compute the interactions at this level directly. */
iact_pair_direct(ci, cj);
}
} else { /* Otherwise, compute the interactions at this level directly. */
iact_pair_direct(ci, cj);
}
}
......@@ -928,7 +941,7 @@ void create_tasks(struct qsched *s, struct cell *ci, struct cell *cj) {
} else {/* Cells are direct neighbours */
/* Are both cells split ? */
/* Are both cells split and we are above the task limit ? */
if (ci->split && cj->split && ci->count > task_limit / cj->count) {
/* Let's split both cells and build all possible pairs */
......@@ -1241,7 +1254,7 @@ void test_bh(int N, int nr_threads, int runs, char *fileName) {
iact_pair(d[0], d[1]);
break;
case task_type_self_pc:
iact_self_pc(d[0], d[1], NULL);
//iact_self_pc(d[0], d[1], NULL);
break;
case task_type_com:
comp_com(d[0]);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment