diff --git a/examples/test_bh_mpi.c b/examples/test_bh_mpi.c
index c6115f32d430300635dec984e8b3091b612e7b50..9bf72a527e7bbbcc79807bdddde2b620e4f09228 100644
--- a/examples/test_bh_mpi.c
+++ b/examples/test_bh_mpi.c
@@ -497,6 +497,11 @@ static inline void iact_pair_pc(struct qsched *s, struct cell *ci, struct cell *
int part_count;
ci->parts = (struct part*) qsched_getresdata(s, ci->res_parts);
+ #ifdef SANITY_CHECKS
+ if(are_neighbours(ci, cj))
+ error("ci and cj are neighbours");
+ #endif
+
/* Load particles. */
part_count = ci->count;
if ((parts = (struct part_local *) malloc(sizeof(struct part_local) * part_count)) == NULL)
@@ -535,104 +540,6 @@ static inline void iact_pair_pc(struct qsched *s, struct cell *ci, struct cell *
free(parts);
}
-#ifdef OLD_VERSION
-/**
- * @brief Interacts a cell ci with the monopole of cell cj, or recurses cj to the depth of ci if ci and cj are neighbours.
- *
- * @param ci The first cell (particles).
- * @param cj The second cell (monopole).
-*/
-static inline void iact_pair_pc(struct qsched *s, struct cell *ci, struct cell *cj) {
-
- struct part_local *parts = NULL;
- struct multipole_local mult_local[8];
- struct cell *cps;
- int part_count, mult_count;
- ci->parts = (struct part* ) qsched_getresdata(s, ci->res_parts);
-
- /* Check if ci and cj are neighbours. */
- if(are_neighbours_different_size(ci, cj) && ci->h != cj->h && cj->split)
- {
- /* They are neighbours, so we need to recurse on cj.*/
- for (cps = (struct cell*) qsched_getresdata(s, cj->firstchild); cps != (struct cell*) qsched_getresdata(s, cj->sibling); cps =
- (struct cell*) qsched_getresdata(s, cps->sibling) ) {
- iact_pair_pc(s, ci, cps);
- }
- }else{ /* Else cj->h >= ci->h and they're not neighbours, so multipole interaction!*/
- /* Load particles. */
- part_count = ci->count;
- if ((parts = (struct part_local *) malloc(sizeof(struct part_local) * part_count)) == NULL)
- error("Failed to allocate local parts.");
- for (int k = 0; k < part_count; k++) {
- for (int j = 0; j < 3; j++) {
- parts[k].x[j] = ci->parts[k].x[j] - ci->loc[j];
- parts[k].a[j] = 0.0f;
- }
- parts[k].mass = ci->parts[k].mass;
- }
- mult_count = 0;
- if(cj->split && cj->h > ci->h)
- {
- /* Now loop over all the other cells */
- for (cps = (struct cell*) qsched_getresdata(s, cj->firstchild); cps != (struct cell*) qsched_getresdata(s, cj->sibling); cps =
- (struct cell*) qsched_getresdata(s, cps->sibling) ) {
-
- mult_local[mult_count].com[0] = cps->new.com[0] - ci->loc[0];
- mult_local[mult_count].com[1] = cps->new.com[1] - ci->loc[1];
- mult_local[mult_count].com[2] = cps->new.com[2] - ci->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_xx = cps->new.I_xx*imass - cps->new.com[0] * cps->new.com[0];
- mult_local[mult_count].I_xy = cps->new.I_xy*imass - cps->new.com[0] * cps->new.com[1];
- mult_local[mult_count].I_xz = cps->new.I_xz*imass - cps->new.com[0] * cps->new.com[2];
- mult_local[mult_count].I_yy = cps->new.I_yy*imass - cps->new.com[1] * cps->new.com[1];
- mult_local[mult_count].I_yz = cps->new.I_yz*imass - cps->new.com[1] * cps->new.com[2];
- mult_local[mult_count].I_zz = cps->new.I_zz*imass - cps->new.com[2] * cps->new.com[2];
- #endif
-
- mult_count++;
-
- }
- }else
- {
- mult_count = 1;
- mult_local[mult_count].com[0] = cj->new.com[0] - ci->loc[0];
- mult_local[mult_count].com[1] = cj->new.com[1] - ci->loc[1];
- mult_local[mult_count].com[2] = cj->new.com[2] - ci->loc[2];
- mult_local[mult_count].mass = cj->new.mass;
-
- #ifdef QUADRUPOLES
-
- /* Final operation to get the quadrupoles */
- double imass = 1. / cj->new.mass;
- mult_local[mult_count].I_xx = cj->new.I_xx*imass - cj->new.com[0] * cj->new.com[0];
- mult_local[mult_count].I_xy = cj->new.I_xy*imass - cj->new.com[0] * cj->new.com[1];
- mult_local[mult_count].I_xz = cj->new.I_xz*imass - cj->new.com[0] * cj->new.com[2];
- mult_local[mult_count].I_yy = cj->new.I_yy*imass - cj->new.com[1] * cj->new.com[1];
- mult_local[mult_count].I_yz = cj->new.I_yz*imass - cj->new.com[1] * cj->new.com[2];
- mult_local[mult_count].I_zz = cj->new.I_zz*imass - cj->new.com[2] * cj->new.com[2];
- #endif
- }
- /* Perform the interaction between the particles and the list of multipoles */
- make_interact_pc(parts, part_count, mult_local, mult_count);
-
- /* Clean up local parts? */
- for (int k = 0; k < part_count; k++) {
- for (int j = 0; j < 3; j++) ci->parts[k].a[j] += parts[k].a[j];
- }
- free(parts);
- }
-
-}
-#endif
-
-
-
-
struct cell *cell_get(struct qsched *s) {
struct cell *res;
@@ -904,6 +811,10 @@ void create_pcs(struct qsched *s, struct cell *ci, struct cell *cj, int depth, i
qsched_task_t cp, cps;
struct cell *cp1, *cp2;
+ #ifdef SANITY_CHECKS
+ if(cj!= NULL && ci->h != cj->h)
+ error("ci and cj have different h values.");
+ #endif
/* If cj == NULL we're at the root.*/
if(cj == NULL)
@@ -1388,6 +1299,7 @@ if(s.rank == 0)
c = (struct cell*) qsched_getresdata( &s, c->firstchild );
}
}
+ message("There are %i leaves.", nr_leaves);
message("Average number of parts per leaf is %f.", ((double)N) / nr_leaves);
}
@@ -1431,7 +1343,9 @@ if(s.rank == 0)
}
}
message("leaf_depth = %i", leaf_depth);
- create_pcs(&s, root, NULL, 0, 2);
+ message("tasks before = %i", s.count);
+ create_pcs(&s, root, NULL, 0, leaf_depth-1);
+ message("tasks after = %i", s.count);
}
printf("s.count = %i\n", s.count);
//Mark all the schedulers dirty so they actually synchronize.