diff --git a/examples/test_fmm.c b/examples/test_fmm.c
index eaa89ee5a0d91c6b85eed1ae39924e843dd49044..33d7d88409ff9492a96afce2295eefc3511f194d 100644
--- a/examples/test_fmm.c
+++ b/examples/test_fmm.c
@@ -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]);