diff --git a/src/cell.h b/src/cell.h
index 747b79f05a799b9401c67c23b72b722d385c5ced..b282a094da7cd53de931638be18acc3e31d565b1 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -282,6 +282,15 @@ struct cell {
/*! Task propagating the multipole to the particles */
struct task *grav_down;
+ /*! Dependency implicit task for the star ghost (in->ghost->out)*/
+ struct task *star_ghost_in;
+
+ /*! Dependency implicit task for the star ghost (in->ghost->out)*/
+ struct task *star_ghost_out;
+
+ /*! The star ghost task itself */
+ struct task *star_ghost;
+
/*! Task for cooling */
struct task *cooling;
diff --git a/src/engine.c b/src/engine.c
index adb853ac91d89c48855ba2b46fc058a5a7805d40..91263f049274070479ad9671571f0a2744d24bba 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -83,6 +83,7 @@
#include "sort_part.h"
#include "sourceterms.h"
#include "statistics.h"
+#include "stars_io.h"
#include "timers.h"
#include "tools.h"
#include "units.h"
@@ -148,6 +149,29 @@ void engine_addlink(struct engine *e, struct link **l, struct task *t) {
res->next = atomic_swap(l, res);
}
+/**
+ * @brief Recursively add non-implicit star ghost tasks to a cell hierarchy.
+ */
+void engine_add_star_ghosts(struct engine *e, struct cell *c, struct task *star_ghost_in,
+ struct task *star_ghost_out) {
+
+ /* If we have reached the leaf OR have to few particles to play with*/
+ if (!c->split || c->scount < engine_max_sparts_per_ghost) {
+
+ /* Add the ghost task and its dependencies */
+ struct scheduler *s = &e->sched;
+ c->star_ghost =
+ scheduler_addtask(s, task_type_star_ghost, task_subtype_none, 0, 0, c, NULL);
+ scheduler_addunlock(s, star_ghost_in, c->star_ghost);
+ scheduler_addunlock(s, c->star_ghost, star_ghost_out);
+ } else {
+ /* Keep recursing */
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL)
+ engine_add_star_ghosts(e, c->progeny[k], star_ghost_in, star_ghost_out);
+ }
+}
+
/**
* @brief Recursively add non-implicit ghost tasks to a cell hierarchy.
*/
@@ -389,6 +413,48 @@ void engine_make_hierarchical_tasks_gravity(struct engine *e, struct cell *c) {
engine_make_hierarchical_tasks_gravity(e, c->progeny[k]);
}
+/**
+ * @brief Generate the stars hierarchical tasks for a hierarchy of cells -
+ * i.e. all the O(Npart) tasks -- star version
+ *
+ * Tasks are only created here. The dependencies will be added later on.
+ *
+ * Note that there is no need to recurse below the super-cell. Note also
+ * that we only add tasks if the relevant particles are present in the cell.
+ *
+ * @param e The #engine.
+ * @param c The #cell.
+ */
+void engine_make_hierarchical_tasks_stars(struct engine *e, struct cell *c) {
+
+ struct scheduler *s = &e->sched;
+
+ /* Are we in a super-cell ? */
+ if (c->super_hydro == c) {
+
+ /* Local tasks only... */
+ if (c->nodeID == e->nodeID) {
+
+ /* Generate the ghost tasks. */
+ c->star_ghost_in =
+ scheduler_addtask(s, task_type_star_ghost_in, task_subtype_none, 0,
+ /* implicit = */ 1, c, NULL);
+ c->star_ghost_out =
+ scheduler_addtask(s, task_type_star_ghost_out, task_subtype_none, 0,
+ /* implicit = */ 1, c, NULL);
+ engine_add_star_ghosts(e, c, c->star_ghost_in, c->star_ghost_out);
+
+ }
+ } else { /* We are above the super-cell so need to go deeper */
+
+ /* Recurse. */
+ if (c->split)
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL)
+ engine_make_hierarchical_tasks_stars(e, c->progeny[k]);
+ }
+}
+
void engine_make_hierarchical_tasks_mapper(void *map_data, int num_elements,
void *extra_data) {
struct engine *e = (struct engine *)extra_data;
@@ -396,6 +462,7 @@ void engine_make_hierarchical_tasks_mapper(void *map_data, int num_elements,
const int is_with_self_gravity = (e->policy & engine_policy_self_gravity);
const int is_with_external_gravity =
(e->policy & engine_policy_external_gravity);
+ const int is_with_stars = (e->policy & engine_policy_stars);
for (int ind = 0; ind < num_elements; ind++) {
struct cell *c = &((struct cell *)map_data)[ind];
@@ -406,6 +473,8 @@ void engine_make_hierarchical_tasks_mapper(void *map_data, int num_elements,
/* And the gravity stuff */
if (is_with_self_gravity || is_with_external_gravity)
engine_make_hierarchical_tasks_gravity(e, c);
+ if (is_with_stars)
+ engine_make_hierarchical_tasks_stars(e, c);
}
}
@@ -2876,7 +2945,22 @@ static inline void engine_make_hydro_loops_dependencies(struct scheduler *sched,
scheduler_addunlock(sched, c->super_hydro->ghost_out, force);
}
+
#endif
+/**
+ * @brief Creates the dependency network for the stars tasks of a given cell.
+ *
+ * @param sched The #scheduler.
+ * @param density The density task to link.
+ * @param c The cell.
+ */
+static inline void engine_make_stars_loops_dependencies(struct scheduler *sched,
+ struct task *density,
+ struct cell *c) {
+ /* density loop --> ghost */
+ scheduler_addunlock(sched, density, c->super_hydro->star_ghost_in);
+}
+
/**
* @brief Duplicates the first hydro loop and construct all the
* dependencies for the hydro part
@@ -2887,13 +2971,14 @@ static inline void engine_make_hydro_loops_dependencies(struct scheduler *sched,
* With all the relevant tasks for a given cell available, we construct
* all the dependencies for that cell.
*/
-void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
+void engine_make_extra_loop_tasks_mapper(void *map_data, int num_elements,
void *extra_data) {
struct engine *e = (struct engine *)extra_data;
struct scheduler *sched = &e->sched;
const int nodeID = e->nodeID;
const int with_cooling = (e->policy & engine_policy_cooling);
+ const int with_stars = (e->policy & engine_policy_stars);
for (int ind = 0; ind < num_elements; ind++) {
struct task *t = &((struct task *)map_data)[ind];
@@ -2923,6 +3008,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* Now, build all the dependencies for the hydro */
engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->ci,
with_cooling);
+ if (with_stars)
+ engine_make_stars_loops_dependencies(sched, t, t->ci);
scheduler_addunlock(sched, t3, t->ci->super->end_force);
#else
@@ -2934,7 +3021,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
engine_addlink(e, &t->ci->force, t2);
/* Now, build all the dependencies for the hydro */
- engine_make_hydro_loops_dependencies(sched, t, t2, t->ci, with_cooling);
+ if (with_stars)
+ engine_make_stars_loops_dependencies(sched, t, t->ci);
scheduler_addunlock(sched, t2, t->ci->super->end_force);
#endif
}
@@ -2970,12 +3058,18 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
if (t->ci->nodeID == nodeID) {
engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->ci,
with_cooling);
+ if (with_stars)
+ engine_make_stars_loops_dependencies(sched, t, t->ci);
+
scheduler_addunlock(sched, t3, t->ci->super->end_force);
}
if (t->cj->nodeID == nodeID) {
if (t->ci->super_hydro != t->cj->super_hydro)
engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->cj,
with_cooling);
+ if (with_stars)
+ engine_make_stars_loops_dependencies(sched, t, t->ci);
+
if (t->ci->super != t->cj->super)
scheduler_addunlock(sched, t3, t->cj->super->end_force);
}
@@ -2994,12 +3088,16 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* that are local and are not descendant of the same super_hydro-cells */
if (t->ci->nodeID == nodeID) {
engine_make_hydro_loops_dependencies(sched, t, t2, t->ci, with_cooling);
+ if (with_stars)
+ engine_make_stars_loops_dependencies(sched, t, t->ci);
scheduler_addunlock(sched, t2, t->ci->super->end_force);
}
if (t->cj->nodeID == nodeID) {
if (t->ci->super_hydro != t->cj->super_hydro)
engine_make_hydro_loops_dependencies(sched, t, t2, t->cj,
with_cooling);
+ if (with_stars)
+ engine_make_stars_loops_dependencies(sched, t, t->ci);
if (t->ci->super != t->cj->super)
scheduler_addunlock(sched, t2, t->cj->super->end_force);
}
@@ -3035,6 +3133,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
if (t->ci->nodeID == nodeID) {
engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->ci,
with_cooling);
+ if (with_stars)
+ engine_make_stars_loops_dependencies(sched, t, t->ci);
scheduler_addunlock(sched, t3, t->ci->super->end_force);
}
@@ -3051,6 +3151,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* that are local and are not descendant of the same super_hydro-cells */
if (t->ci->nodeID == nodeID) {
engine_make_hydro_loops_dependencies(sched, t, t2, t->ci, with_cooling);
+ if (with_stars)
+ engine_make_stars_loops_dependencies(sched, t, t->ci);
scheduler_addunlock(sched, t2, t->ci->super->end_force);
}
#endif
@@ -3091,12 +3193,16 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
if (t->ci->nodeID == nodeID) {
engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->ci,
with_cooling);
+ if (with_stars)
+ engine_make_stars_loops_dependencies(sched, t, t->ci);
scheduler_addunlock(sched, t3, t->ci->super->end_force);
}
if (t->cj->nodeID == nodeID) {
if (t->ci->super_hydro != t->cj->super_hydro)
engine_make_hydro_loops_dependencies(sched, t, t2, t3, t->cj,
with_cooling);
+ if (with_stars)
+ engine_make_stars_loops_dependencies(sched, t, t->ci);
if (t->ci->super != t->cj->super)
scheduler_addunlock(sched, t3, t->cj->super->end_force);
}
@@ -3115,12 +3221,16 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* that are local and are not descendant of the same super_hydro-cells */
if (t->ci->nodeID == nodeID) {
engine_make_hydro_loops_dependencies(sched, t, t2, t->ci, with_cooling);
+ if (with_stars)
+ engine_make_stars_loops_dependencies(sched, t, t->ci);
scheduler_addunlock(sched, t2, t->ci->super->end_force);
}
if (t->cj->nodeID == nodeID) {
if (t->ci->super_hydro != t->cj->super_hydro)
engine_make_hydro_loops_dependencies(sched, t, t2, t->cj,
with_cooling);
+ if (with_stars)
+ engine_make_stars_loops_dependencies(sched, t, t->ci);
if (t->ci->super != t->cj->super)
scheduler_addunlock(sched, t2, t->cj->super->end_force);
}
diff --git a/src/engine.h b/src/engine.h
index 5bb2780211f953edefa192e90e4ae1076409c4a5..3b09abcf112ccbaaf53ee5d4ac192c1c983481b0 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -97,6 +97,7 @@ enum engine_step_properties {
#define engine_default_energy_file_name "energy"
#define engine_default_timesteps_file_name "timesteps"
#define engine_max_parts_per_ghost 1000
+#define engine_max_sparts_per_ghost 1000
/**
* @brief The rank of the engine as a global variable (for messages).
@@ -394,7 +395,8 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
const struct unit_system *internal_units,
const struct phys_const *physical_constants,
struct cosmology *cosmo, const struct hydro_props *hydro,
- struct gravity_props *gravity, struct pm_mesh *mesh,
+ struct gravity_props *gravity, const struct stars_props *stars,
+ struct pm_mesh *mesh,
const struct external_potential *potential,
const struct cooling_function_data *cooling_func,
const struct chemistry_global_data *chemistry,
diff --git a/src/runner.c b/src/runner.c
index cbd4e6d7535037356997b286d256b740697388f9..5692e7d1e433fabb466a7925ffb98b396aaec20c 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -147,7 +147,6 @@ void runner_do_star_ghost(struct runner *r, struct cell *c, int timer) {
struct spart *restrict sparts = c->sparts;
const struct engine *e = r->e;
- const struct space *s = e->s;
const struct cosmology *cosmo = e->cosmology;
const struct stars_props *stars_properties = e->stars_properties;
const float star_h_max = e->stars_properties->h_max;
@@ -171,7 +170,7 @@ void runner_do_star_ghost(struct runner *r, struct cell *c, int timer) {
/* Init the list of active particles that have to be updated. */
int *sid = NULL;
if ((sid = (int *)malloc(sizeof(int) * c->scount)) == NULL)
- error("Can't allocate memory for pid.");
+ error("Can't allocate memory for sid.");
for (int k = 0; k < c->scount; k++)
if (spart_is_active(&sparts[k], e)) {
sid[count] = k;
@@ -303,23 +302,23 @@ void runner_do_star_ghost(struct runner *r, struct cell *c, int timer) {
/* Self-interaction? */
if (l->t->type == task_type_self)
- runner_doself_subset_star_density(r, finger, parts, pid, count);
+ runner_doself_subset_branch_star_density(r, finger, sparts, sid, count);
/* Otherwise, pair interaction? */
else if (l->t->type == task_type_pair) {
/* Left or right? */
if (l->t->ci == finger)
- runner_dopair_subset_star_density(r, finger, parts, pid,
+ runner_dopair_subset_branch_star_density(r, finger, sparts, sid,
count, l->t->cj);
else
- runner_dopair_subset_star_density(r, finger, parts, pid,
+ runner_dopair_subset_branch_star_density(r, finger, sparts, sid,
count, l->t->ci);
}
/* Otherwise, sub-self interaction? */
else if (l->t->type == task_type_sub_self)
- runner_dosub_subset_density(r, finger, parts, pid, count, NULL,
+ runner_dosub_subset_star_density(r, finger, sparts, sid, count, NULL,
-1, 1);
/* Otherwise, sub-pair interaction? */
@@ -327,10 +326,10 @@ void runner_do_star_ghost(struct runner *r, struct cell *c, int timer) {
/* Left or right? */
if (l->t->ci == finger)
- runner_dosub_subset_density(r, finger, parts, pid, count,
+ runner_dosub_subset_star_density(r, finger, sparts, sid, count,
l->t->cj, -1, 1);
else
- runner_dosub_subset_density(r, finger, parts, pid, count,
+ runner_dosub_subset_star_density(r, finger, sparts, sid, count,
l->t->ci, -1, 1);
}
}
diff --git a/src/runner_doiact_star.h b/src/runner_doiact_star.h
index 555b665c81b88380c9c1a1a96d292c521a203ed1..e02b8f157dd0234e5b77b66916a00adf145ce74b 100644
--- a/src/runner_doiact_star.h
+++ b/src/runner_doiact_star.h
@@ -251,3 +251,698 @@ void runner_dopair_subset_star_density(struct runner *r, struct cell *restrict c
TIMER_TOC(timer_dopair_subset_naive);
}
+
+
+/**
+ * @brief Compute the interactions between a cell pair, but only for the
+ * given indices in ci.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param sparts The #spart to interact.
+ * @param ind The list of indices of particles in @c ci to interact with.
+ * @param scount The number of particles in @c ind.
+ */
+void runner_doself_subset_star_density(struct runner *r, struct cell *restrict ci,
+ struct spart *restrict sparts, int *restrict ind, int scount) {
+
+ const struct engine *e = r->e;
+ const struct cosmology *cosmo = e->cosmology;
+
+ TIMER_TIC;
+
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
+ const int count_i = ci->count;
+ struct part *restrict parts_j = ci->parts;
+
+ /* Loop over the parts in ci. */
+ for (int spid = 0; spid < scount; spid++) {
+
+ /* Get a hold of the ith part in ci. */
+ struct spart *spi = &sparts[ind[spid]];
+ const float spix[3] = {(float)(spi->x[0] - ci->loc[0]),
+ (float)(spi->x[1] - ci->loc[1]),
+ (float)(spi->x[2] - ci->loc[2])};
+ const float hi = spi->h;
+ const float hig2 = hi * hi * kernel_gamma2;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (!spart_is_active(spi, e)) error("Inactive particle in subset function!");
+#endif
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count_i; pjd++) {
+
+ /* Get a pointer to the jth particle. */
+ struct part *restrict pj = &parts_j[pjd];
+ const float hj = pj->h;
+
+ /* Compute the pairwise distance. */
+ const float pjx[3] = {(float)(pj->x[0] - ci->loc[0]),
+ (float)(pj->x[1] - ci->loc[1]),
+ (float)(pj->x[2] - ci->loc[2])};
+ float dx[3] = {spix[0] - pjx[0], spix[1] - pjx[1], spix[2] - pjx[2]};
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles have been drifted to the current time */
+ if (spi->ti_drift != e->ti_current)
+ error("Particle pi not drifted to current time");
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
+
+ /* Hit or miss? */
+ if (r2 > 0.f && r2 < hig2) {
+ runner_iact_nonsym_star_density(r2, dx, hi, hj, spi, pj, a, H);
+ }
+ } /* loop over the parts in cj. */
+ } /* loop over the parts in ci. */
+
+ TIMER_TOC(timer_doself_subset_star_density);
+}
+
+
+ /**
+ * @brief Determine which version of DOSELF_SUBSET needs to be called depending
+ * on the optimisation level.
+
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param parts The #spart to interact.
+ * @param ind The list of indices of particles in @c ci to interact with.
+ * @param scount The number of particles in @c ind.
+ */
+void runner_doself_subset_branch_star_density(struct runner *r, struct cell *restrict ci,
+ struct spart *restrict sparts, int *restrict ind,
+ int scount) {
+
+ runner_doself_subset_star_density(r, ci, sparts, ind, scount);
+}
+
+ /**
+ * @brief Determine which version of DOPAIR_SUBSET needs to be called depending
+ * on the
+ * orientation of the cells or whether DOPAIR_SUBSET needs to be called at all.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param sparts_i The #spart to interact with @c cj.
+ * @param ind The list of indices of particles in @c ci to interact with.
+ * @param scount The number of particles in @c ind.
+ * @param cj The second #cell.
+ */
+ void runner_dopair_subset_branch_star_density(struct runner *r, struct cell *restrict ci,
+ struct spart *restrict sparts_i, int *restrict ind,
+ int scount, struct cell *restrict cj) {
+
+ const struct engine *e = r->e;
+
+ /* Get the relative distance between the pairs, wrapping. */
+ double shift[3] = {0.0, 0.0, 0.0};
+ for (int k = 0; k < 3; k++) {
+ if (cj->loc[k] - ci->loc[k] < -e->s->dim[k] / 2)
+ shift[k] = e->s->dim[k];
+ else if (cj->loc[k] - ci->loc[k] > e->s->dim[k] / 2)
+ shift[k] = -e->s->dim[k];
+ }
+
+ runner_dopair_subset_star_density(r, ci, sparts_i, ind, scount, cj, shift);
+}
+
+void runner_dosub_subset_star_density(struct runner *r, struct cell *ci, struct spart *sparts,
+ int *ind, int scount, struct cell *cj, int sid, int gettimer) {
+
+ const struct engine *e = r->e;
+ struct space *s = e->s;
+
+ TIMER_TIC;
+
+ /* Should we even bother? */
+ if (!cell_is_active_star(ci, e) &&
+ (cj == NULL || !cell_is_active_star(cj, e)))
+ return;
+ if (ci->scount == 0 || (cj != NULL && cj->scount == 0)) return;
+
+ /* Find out in which sub-cell of ci the parts are. */
+ struct cell *sub = NULL;
+ if (ci->split) {
+ for (int k = 0; k < 8; k++) {
+ if (ci->progeny[k] != NULL) {
+ if (&sparts[ind[0]] >= &ci->progeny[k]->sparts[0] &&
+ &sparts[ind[0]] < &ci->progeny[k]->sparts[ci->progeny[k]->scount]) {
+ sub = ci->progeny[k];
+ break;
+ }
+ }
+ }
+ }
+
+ /* Is this a single cell? */
+ if (cj == NULL) {
+
+ /* Recurse? */
+ if (cell_can_recurse_in_self_task(ci)) {
+
+ /* Loop over all progeny. */
+ runner_dosub_subset_star_density(r, sub, sparts, ind, scount, NULL, -1, 0);
+ for (int j = 0; j < 8; j++)
+ if (ci->progeny[j] != sub && ci->progeny[j] != NULL)
+ runner_dosub_subset_star_density(r, sub, sparts, ind, scount, ci->progeny[j], -1, 0);
+
+ }
+
+ /* Otherwise, compute self-interaction. */
+ else
+ runner_doself_subset_branch_star_density(r, ci, sparts, ind, scount);
+ } /* self-interaction. */
+
+ /* Otherwise, it's a pair interaction. */
+ else {
+
+ /* Recurse? */
+ if (cell_can_recurse_in_pair_task(ci) &&
+ cell_can_recurse_in_pair_task(cj)) {
+
+ /* Get the type of pair if not specified explicitly. */
+ double shift[3] = {0.0, 0.0, 0.0};
+ sid = space_getsid(s, &ci, &cj, shift);
+
+ /* Different types of flags. */
+ switch (sid) {
+
+ /* Regular sub-cell interactions of a single cell. */
+ case 0: /* ( 1 , 1 , 1 ) */
+ if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ break;
+
+ case 1: /* ( 1 , 1 , 0 ) */
+ if (ci->progeny[6] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ break;
+
+ case 2: /* ( 1 , 1 , -1 ) */
+ if (ci->progeny[6] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ break;
+
+ case 3: /* ( 1 , 0 , 1 ) */
+ if (ci->progeny[5] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ if (ci->progeny[5] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ break;
+
+ case 4: /* ( 1 , 0 , 0 ) */
+ if (ci->progeny[4] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[4] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[4],
+ -1, 0);
+ if (ci->progeny[4] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[4] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[4],
+ -1, 0);
+ if (ci->progeny[4] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[4] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[4],
+ -1, 0);
+ if (ci->progeny[4] == sub && cj->progeny[3] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[3],
+ -1, 0);
+ if (ci->progeny[4] != NULL && cj->progeny[3] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[4],
+ -1, 0);
+ if (ci->progeny[5] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ if (ci->progeny[5] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ if (ci->progeny[5] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ if (ci->progeny[5] == sub && cj->progeny[3] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[3],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[3] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[3] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[3],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[3] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[3] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[3],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[3] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ break;
+
+ case 5: /* ( 1 , 0 , -1 ) */
+ if (ci->progeny[4] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[4] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[4],
+ -1, 0);
+ if (ci->progeny[4] == sub && cj->progeny[3] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[3],
+ -1, 0);
+ if (ci->progeny[4] != NULL && cj->progeny[3] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[4],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[3] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[3],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[3] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ break;
+
+ case 6: /* ( 1 , -1 , 1 ) */
+ if (ci->progeny[5] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ break;
+
+ case 7: /* ( 1 , -1 , 0 ) */
+ if (ci->progeny[4] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[4] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[4],
+ -1, 0);
+ if (ci->progeny[4] == sub && cj->progeny[3] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[3],
+ -1, 0);
+ if (ci->progeny[4] != NULL && cj->progeny[3] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[4],
+ -1, 0);
+ if (ci->progeny[5] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ if (ci->progeny[5] == sub && cj->progeny[3] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[3],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[3] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ break;
+
+ case 8: /* ( 1 , -1 , -1 ) */
+ if (ci->progeny[4] == sub && cj->progeny[3] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[4], sparts, ind, scount, cj->progeny[3],
+ -1, 0);
+ if (ci->progeny[4] != NULL && cj->progeny[3] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[3], sparts, ind, scount, ci->progeny[4],
+ -1, 0);
+ break;
+
+ case 9: /* ( 0 , 1 , 1 ) */
+ if (ci->progeny[3] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[3] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[3],
+ -1, 0);
+ if (ci->progeny[3] == sub && cj->progeny[4] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[4],
+ -1, 0);
+ if (ci->progeny[3] != NULL && cj->progeny[4] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[3],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[4] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[4],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[4] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ break;
+
+ case 10: /* ( 0 , 1 , 0 ) */
+ if (ci->progeny[2] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[2] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[2],
+ -1, 0);
+ if (ci->progeny[2] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[2] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[2],
+ -1, 0);
+ if (ci->progeny[2] == sub && cj->progeny[4] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[4],
+ -1, 0);
+ if (ci->progeny[2] != NULL && cj->progeny[4] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[2],
+ -1, 0);
+ if (ci->progeny[2] == sub && cj->progeny[5] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[5],
+ -1, 0);
+ if (ci->progeny[2] != NULL && cj->progeny[5] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[2],
+ -1, 0);
+ if (ci->progeny[3] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[3] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[3],
+ -1, 0);
+ if (ci->progeny[3] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[3] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[3],
+ -1, 0);
+ if (ci->progeny[3] == sub && cj->progeny[4] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[4],
+ -1, 0);
+ if (ci->progeny[3] != NULL && cj->progeny[4] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[3],
+ -1, 0);
+ if (ci->progeny[3] == sub && cj->progeny[5] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[5],
+ -1, 0);
+ if (ci->progeny[3] != NULL && cj->progeny[5] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[3],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[4] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[4],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[4] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[5] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[5],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[5] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[4] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[4],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[4] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[5] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[5],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[5] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ break;
+
+ case 11: /* ( 0 , 1 , -1 ) */
+ if (ci->progeny[2] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[2] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[2],
+ -1, 0);
+ if (ci->progeny[2] == sub && cj->progeny[5] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[2], sparts, ind, scount, cj->progeny[5],
+ -1, 0);
+ if (ci->progeny[2] != NULL && cj->progeny[5] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[2],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[1] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[1],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[1] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[1], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ if (ci->progeny[6] == sub && cj->progeny[5] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[6], sparts, ind, scount, cj->progeny[5],
+ -1, 0);
+ if (ci->progeny[6] != NULL && cj->progeny[5] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[5], sparts, ind, scount, ci->progeny[6],
+ -1, 0);
+ break;
+
+ case 12: /* ( 0 , 0 , 1 ) */
+ if (ci->progeny[1] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[1], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[1] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[1],
+ -1, 0);
+ if (ci->progeny[1] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[1], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[1] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[1],
+ -1, 0);
+ if (ci->progeny[1] == sub && cj->progeny[4] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[1], sparts, ind, scount, cj->progeny[4],
+ -1, 0);
+ if (ci->progeny[1] != NULL && cj->progeny[4] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[1],
+ -1, 0);
+ if (ci->progeny[1] == sub && cj->progeny[6] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[1], sparts, ind, scount, cj->progeny[6],
+ -1, 0);
+ if (ci->progeny[1] != NULL && cj->progeny[6] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[6], sparts, ind, scount, ci->progeny[1],
+ -1, 0);
+ if (ci->progeny[3] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[3] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[3],
+ -1, 0);
+ if (ci->progeny[3] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[3] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[3],
+ -1, 0);
+ if (ci->progeny[3] == sub && cj->progeny[4] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[4],
+ -1, 0);
+ if (ci->progeny[3] != NULL && cj->progeny[4] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[3],
+ -1, 0);
+ if (ci->progeny[3] == sub && cj->progeny[6] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[3], sparts, ind, scount, cj->progeny[6],
+ -1, 0);
+ if (ci->progeny[3] != NULL && cj->progeny[6] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[6], sparts, ind, scount, ci->progeny[3],
+ -1, 0);
+ if (ci->progeny[5] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ if (ci->progeny[5] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ if (ci->progeny[5] == sub && cj->progeny[4] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[4],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[4] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ if (ci->progeny[5] == sub && cj->progeny[6] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[5], sparts, ind, scount, cj->progeny[6],
+ -1, 0);
+ if (ci->progeny[5] != NULL && cj->progeny[6] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[6], sparts, ind, scount, ci->progeny[5],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[0] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[0],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[0] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[0], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[2] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[2],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[2] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[2], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[4] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[4],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[4] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[4], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ if (ci->progeny[7] == sub && cj->progeny[6] != NULL)
+ runner_dosub_subset_star_density(r, ci->progeny[7], sparts, ind, scount, cj->progeny[6],
+ -1, 0);
+ if (ci->progeny[7] != NULL && cj->progeny[6] == sub)
+ runner_dosub_subset_star_density(r, cj->progeny[6], sparts, ind, scount, ci->progeny[7],
+ -1, 0);
+ break;
+ }
+
+ }
+
+ /* Otherwise, compute the pair directly. */
+ else if (cell_is_active_hydro(ci, e) || cell_is_active_hydro(cj, e)) {
+
+ /* Do any of the cells need to be drifted first? */
+ if (!cell_are_part_drifted(cj, e)) error("Cell should be drifted!");
+
+ runner_dopair_subset_branch_star_density(r, ci, sparts, ind, scount, cj);
+ }
+
+ } /* otherwise, pair interaction. */
+
+ if (gettimer) TIMER_TOC(timer_dosub_subset);
+}
diff --git a/src/scheduler.c b/src/scheduler.c
index a90a3f1b0fa2401a6bc128069486c679df535435..4c929458e13c2b54e0602f5c108a4edf6fa4d651 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1525,6 +1525,7 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
break;
case task_type_sort:
case task_type_ghost:
+ case task_type_star_ghost:
case task_type_drift_part:
qid = t->ci->super_hydro->owner;
break;
diff --git a/src/stars/Default/star_io.h b/src/stars/Default/star_io.h
index 3c424dc0ae703896605ccb8f23b7b9684f0838b9..6c77315d7b72090a262b3112db95cfaa5a619cfb 100644
--- a/src/stars/Default/star_io.h
+++ b/src/stars/Default/star_io.h
@@ -33,7 +33,7 @@ INLINE static void star_read_particles(struct spart* sparts,
struct io_props* list, int* num_fields) {
/* Say how much we want to read */
- *num_fields = 4;
+ *num_fields = 5;
/* List what we want to read */
list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
@@ -44,6 +44,8 @@ INLINE static void star_read_particles(struct spart* sparts,
sparts, mass);
list[3] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY,
UNIT_CONV_NO_UNITS, sparts, id);
+ list[4] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_LENGTH, sparts, h);
}
/**
@@ -58,7 +60,7 @@ INLINE static void star_write_particles(const struct spart* sparts,
int* num_fields) {
/* Say how much we want to read */
- *num_fields = 4;
+ *num_fields = 5;
/* List what we want to read */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -69,6 +71,8 @@ INLINE static void star_write_particles(const struct spart* sparts,
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, sparts, mass);
list[3] = io_make_output_field("ParticleIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS,
sparts, id);
+ list[4] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+ sparts, h);
}
/**
@@ -122,7 +126,7 @@ INLINE static void stars_props_init(struct stars_props *sp,
*
* @param p The #stars_props.
*/
-void stars_props_print(const struct stars_props *sp) {
+INLINE static void stars_props_print(const struct stars_props *sp) {
/* Now stars */
message("Stars kernel: %s with eta=%f (%.2f neighbours).", kernel_name,
@@ -144,7 +148,7 @@ void stars_props_print(const struct stars_props *sp) {
}
#if defined(HAVE_HDF5)
-void stars_props_print_snapshot(hid_t h_grpstars, const struct stars_props *sp) {
+INLINE static void stars_props_print_snapshot(hid_t h_grpstars, const struct stars_props *sp) {
io_write_attribute_s(h_grpstars, "Kernel function", kernel_name);
io_write_attribute_f(h_grpstars, "Kernel target N_ngb", sp->target_neighbours);
@@ -167,7 +171,7 @@ void stars_props_print_snapshot(hid_t h_grpstars, const struct stars_props *sp)
* @param p the struct
* @param stream the file stream
*/
-void stars_props_struct_dump(const struct stars_props *p, FILE *stream) {
+INLINE static void stars_props_struct_dump(const struct stars_props *p, FILE *stream) {
restart_write_blocks((void *)p, sizeof(struct stars_props), 1, stream,
"starsprops", "stars props");
}
@@ -179,7 +183,7 @@ void stars_props_struct_dump(const struct stars_props *p, FILE *stream) {
* @param p the struct
* @param stream the file stream
*/
-void stars_props_struct_restore(const struct stars_props *p, FILE *stream) {
+INLINE static void stars_props_struct_restore(const struct stars_props *p, FILE *stream) {
restart_read_blocks((void *)p, sizeof(struct stars_props), 1, stream, NULL,
"stars props");
}
diff --git a/src/swift.h b/src/swift.h
index e10938addb99956c202b3e4dd2b0592b580fa948..4398f3f69a66320e13d0ed1a6a0a63fafc7d1e52 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -66,6 +66,8 @@
#include "single_io.h"
#include "sourceterms.h"
#include "space.h"
+#include "stars.h"
+#include "stars_io.h"
#include "task.h"
#include "threadpool.h"
#include "timeline.h"
diff --git a/src/task.c b/src/task.c
index 10f2ddf5cec885ec23c4f65db5cdea50f0e5097b..4eb2b0ca7a3abfaf7cf6d9615a294fef17210c89 100644
--- a/src/task.c
+++ b/src/task.c
@@ -48,15 +48,16 @@
/* Task type names. */
const char *taskID_names[task_type_count] = {
- "none", "sort", "self",
- "pair", "sub_self", "sub_pair",
- "init_grav", "init_grav_out", "ghost_in",
- "ghost", "ghost_out", "extra_ghost",
- "drift_part", "drift_gpart", "end_force",
- "kick1", "kick2", "timestep",
- "send", "recv", "grav_long_range",
- "grav_mm", "grav_down_in", "grav_down",
- "grav_mesh", "cooling", "sourceterms"};
+ "none", "sort", "self",
+ "pair", "sub_self", "sub_pair",
+ "init_grav", "init_grav_out", "ghost_in",
+ "ghost", "ghost_out", "extra_ghost",
+ "drift_part", "drift_gpart", "end_force",
+ "kick1", "kick2", "timestep",
+ "send", "recv", "grav_long_range",
+ "grav_mm", "grav_down_in", "grav_down",
+ "grav_mesh", "cooling", "sourceterms",
+ "star_ghost_in", "star_ghost", "star_ghost_out"};
/* Sub-task type names. */
const char *subtaskID_names[task_subtype_count] = {
diff --git a/src/task.h b/src/task.h
index 02c92e1c238c9a92ea1cc6891552420fb3abc0da..eb82b1282960d2f5e058b2940e3b09c19fd6b01b 100644
--- a/src/task.h
+++ b/src/task.h
@@ -66,6 +66,9 @@ enum task_types {
task_type_grav_mesh,
task_type_cooling,
task_type_sourceterms,
+ task_type_star_ghost_in,
+ task_type_star_ghost,
+ task_type_star_ghost_out,
task_type_count
} __attribute__((packed));