Skip to content
Snippets Groups Projects
Commit 0fafeb3f authored by d74ksy's avatar d74ksy
Browse files

Checking this is up to date

parent 28dd9631
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 13 additions and 99 deletions
examples/test_bh_mpi.c
+ 13
− 99
View file @ 0fafeb3f
......@@ -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.
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment