diff --git a/examples/test_bh_new.c b/examples/test_bh_new.c index 2e838379b3e9ffc0a4001c22f49d9339179c57a7..76be2f978e548a9b6246b1d182529fd32b5e276e 100644 --- a/examples/test_bh_new.c +++ b/examples/test_bh_new.c @@ -38,14 +38,14 @@ /* Some local constants. */ #define cell_pool_grow 1000 #define cell_maxparts 100 -#define task_limit 1e8 +#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 NO_COM_AS_TASK +#define COM_AS_TASK #define NO_COUNTERS /** Data structure for the particles. */ @@ -416,7 +416,8 @@ static inline int are_neighbours(struct cell *ci, struct cell *cj) { * @param count The number of particles in the array * @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_local *parts, int count, + double *loc, struct cell *cj) { int j, k; float com[3] = {0.0, 0.0, 0.0}; @@ -441,7 +442,7 @@ static inline void make_interact_pc(struct part *parts, int count, struct cell * #endif /* 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] - loc[k]; mcom = cj->new.mass; /* Loop over every particle in leaf. */ @@ -473,7 +474,9 @@ 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) { + struct part_local *parts = NULL; struct cell *cp, *cps; + int count; #ifdef SANITY_CHECKS if (ci->h != cj->h) @@ -483,17 +486,41 @@ static inline void iact_pair_pc(struct cell *ci, struct cell *cj) { error(" Cells that are not neighbours !"); #endif - /* Let's split both cells and build all possible non-neighbouring pairs */ + /* Let's loop over all siblings of ci */ for (cp = ci->firstchild; cp != ci->sibling; cp = cp->sibling) { + + count = cp->count; + + /* Get local copies of the particle data in cp. */ + 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] = cp->parts[k].x[j] - cp->loc[j]; + parts[k].a[j] = 0.0f; + } + parts[k].mass = cp->parts[k].mass; + } + + + /* Now loop over all the other cells */ for (cps = cj->firstchild; cps != cj->sibling; cps = cps->sibling) { + /* And interact with the ones that are not neighbours */ if ( ! are_neighbours(cp, cps) ) { - make_interact_pc(cp->parts, cp->count, cps); - make_interact_pc(cps->parts, cps->count, cp ); + make_interact_pc(parts, cp->count, cp->loc, cps); } } + + + /* Clean up local parts? */ + for (int k = 0; k < count; k++) { + for (int j = 0; j < 3; j++) cp->parts[k].a[j] += parts[k].a[j]; + } + free(parts); } } @@ -648,6 +675,7 @@ void iact_pair(struct cell *ci, struct cell *cj) { /* Build all possible cell-multipole interaction pairs */ iact_pair_pc(ci, cj); + iact_pair_pc(cj, ci); } else { /* Otherwise, compute the interactions at this level directly. */ iact_pair_direct(ci, cj); @@ -829,23 +857,34 @@ void create_tasks(struct qsched *s, struct cell *ci, struct cell *cj) { } } - /* Let's also build a particle-monopole task */ + /* Let's also build particle-monopole tasks */ /* Create the task. */ data[0] = ci; data[1] = cj; tid = qsched_addtask(s, task_type_pair_pc, task_flag_none, data, - sizeof(struct cell *) * 2, ci->count + cj->count); + sizeof(struct cell *) * 2, ci->count * 8 ); /* Add the resource and dependance */ qsched_addlock(s, tid, ci->res); - qsched_addlock(s, tid, cj->res); #ifdef COM_AS_TASK qsched_addunlock(s, cj->com_tid, tid); - qsched_addunlock(s, ci->com_tid, tid); #endif + + /* Create the task. */ + data[0] = cj; + data[1] = ci; + tid = qsched_addtask(s, task_type_pair_pc, task_flag_none, data, + sizeof(struct cell *) * 2, cj->count * 8); + + /* Add the resource and dependance */ + qsched_addunlock(s, ci->com_tid, tid); + +#ifdef COM_AS_TASK + qsched_addlock(s, tid, cj->res); +#endif } else { /* Otherwise, at least one of the cells is not split, build a direct * interaction. */