From a3505366e5677015c8f08ce35a9eb16fff23c139 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Sat, 16 Feb 2019 21:42:15 +0100
Subject: [PATCH] Added a separate function to drift the stars independantly of
the gparts.
---
src/cell.c | 119 +++++++++++++++++++++++++++++++-------
src/cell.h | 16 ++++-
src/runner.c | 18 ++++++
src/runner_doiact_stars.h | 1 +
4 files changed, 131 insertions(+), 23 deletions(-)
diff --git a/src/cell.c b/src/cell.c
index 22ddc6c520..990823bd9a 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -3830,15 +3830,9 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
const int periodic = e->s->periodic;
const double dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
const int with_cosmology = (e->policy & engine_policy_cosmology);
- const float stars_h_max = e->stars_properties->h_max;
const integertime_t ti_old_gpart = c->grav.ti_old_part;
const integertime_t ti_current = e->ti_current;
struct gpart *const gparts = c->grav.parts;
- struct spart *const sparts = c->stars.parts;
-
- float dx_max = 0.f, dx2_max = 0.f;
- float dx_max_sort = 0.0f, dx2_max_sort = 0.f;
- float cell_h_max = 0.f;
/* Drift irrespective of cell flags? */
force |= c->grav.do_drift;
@@ -3876,19 +3870,9 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
/* Recurse */
cell_drift_gpart(cp, e, force);
-
- /* Update */
- dx_max = max(dx_max, cp->stars.dx_max_part);
- dx_max_sort = max(dx_max_sort, cp->stars.dx_max_sort);
- cell_h_max = max(cell_h_max, cp->stars.h_max);
}
}
- /* Store the values */
- c->stars.h_max = cell_h_max;
- c->stars.dx_max_part = dx_max;
- c->stars.dx_max_sort = dx_max_sort;
-
/* Update the time of the last drift */
c->grav.ti_old_part = ti_current;
@@ -3946,6 +3930,99 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
}
}
+ /* Update the time of the last drift */
+ c->grav.ti_old_part = ti_current;
+ }
+
+ /* Clear the drift flags. */
+ c->grav.do_drift = 0;
+ c->grav.do_sub_drift = 0;
+}
+
+/**
+ * @brief Recursively drifts the #spart in a cell hierarchy.
+ *
+ * @param c The #cell.
+ * @param e The #engine (to get ti_current).
+ * @param force Drift the particles irrespective of the #cell flags.
+ */
+void cell_drift_spart(struct cell *c, const struct engine *e, int force) {
+
+ const int periodic = e->s->periodic;
+ const double dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const float stars_h_max = e->stars_properties->h_max;
+ const integertime_t ti_old_spart = c->stars.ti_old_part;
+ const integertime_t ti_current = e->ti_current;
+ struct spart *const sparts = c->stars.parts;
+
+ float dx_max = 0.f, dx2_max = 0.f;
+ float dx_max_sort = 0.0f, dx2_max_sort = 0.f;
+ float cell_h_max = 0.f;
+
+ /* Drift irrespective of cell flags? */
+ force |= c->stars.do_drift;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that we only drift local cells. */
+ if (c->nodeID != engine_rank) error("Drifting a foreign cell is nope.");
+
+ /* Check that we are actually going to move forward. */
+ if (ti_current < ti_old_spart) error("Attempt to drift to the past");
+#endif
+
+ /* Early abort? */
+ if (c->stars.count == 0) {
+
+ /* Clear the drift flags. */
+ c->stars.do_drift = 0;
+ c->stars.do_sub_drift = 0;
+
+ /* Update the time of the last drift */
+ c->stars.ti_old_part = ti_current;
+
+ return;
+ }
+
+ /* Ok, we have some particles somewhere in the hierarchy to drift */
+
+ /* Are we not in a leaf ? */
+ if (c->split && (force || c->stars.do_sub_drift)) {
+
+ /* Loop over the progeny and collect their data. */
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+ struct cell *cp = c->progeny[k];
+
+ /* Recurse */
+ cell_drift_gpart(cp, e, force);
+
+ /* Update */
+ dx_max = max(dx_max, cp->stars.dx_max_part);
+ dx_max_sort = max(dx_max_sort, cp->stars.dx_max_sort);
+ cell_h_max = max(cell_h_max, cp->stars.h_max);
+ }
+ }
+
+ /* Store the values */
+ c->stars.h_max = cell_h_max;
+ c->stars.dx_max_part = dx_max;
+ c->stars.dx_max_sort = dx_max_sort;
+
+ /* Update the time of the last drift */
+ c->stars.ti_old_part = ti_current;
+
+ } else if (!c->split && force && ti_current > ti_old_spart) {
+
+ /* Drift from the last time the cell was drifted to the current time */
+ double dt_drift;
+ if (with_cosmology) {
+ dt_drift =
+ cosmology_get_drift_factor(e->cosmology, ti_old_spart, ti_current);
+ } else {
+ dt_drift = (ti_current - ti_old_spart) * e->time_base;
+ }
+
/* Loop over all the star particles in the cell */
const size_t nr_sparts = c->stars.count;
for (size_t k = 0; k < nr_sparts; k++) {
@@ -3957,7 +4034,7 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
if (spart_is_inhibited(sp, e)) continue;
/* Drift... */
- drift_spart(sp, dt_drift, ti_old_gpart, ti_current);
+ drift_spart(sp, dt_drift, ti_old_spart, ti_current);
#ifdef SWIFT_DEBUG_CHECKS
/* Make sure the particle does not drift by more than a box length. */
@@ -4009,7 +4086,7 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
if (spart_is_active(sp, e)) {
stars_init_spart(sp);
}
- } /* Note: no need to compute dx_max as all spart have a gpart */
+ }
/* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max);
@@ -4021,12 +4098,12 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
c->stars.dx_max_sort = dx_max_sort;
/* Update the time of the last drift */
- c->grav.ti_old_part = ti_current;
+ c->stars.ti_old_part = ti_current;
}
/* Clear the drift flags. */
- c->grav.do_drift = 0;
- c->grav.do_sub_drift = 0;
+ c->stars.do_drift = 0;
+ c->stars.do_sub_drift = 0;
}
/**
diff --git a/src/cell.h b/src/cell.h
index 63937ba021..a299af5a43 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -490,6 +490,9 @@ struct cell {
/*! Max smoothing length in this cell. */
double h_max;
+ /*! Last (integer) time the cell's spart were drifted forward in time. */
+ integertime_t ti_old_part;
+
/*! Spin lock for various uses (#spart case). */
swift_lock_type lock;
@@ -541,6 +544,12 @@ struct cell {
/*! Is the #spart data of this cell being used in a sub-cell? */
int hold;
+ /*! Does this cell need to be drifted (stars)? */
+ char do_drift;
+
+ /*! Do any of this cell's sub-cells need to be drifted (stars)? */
+ char do_sub_drift;
+
#ifdef SWIFT_DEBUG_CHECKS
/*! Last (integer) time the cell's sort arrays were updated. */
integertime_t ti_sort;
@@ -720,6 +729,7 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s);
void cell_set_super(struct cell *c, struct cell *super);
void cell_drift_part(struct cell *c, const struct engine *e, int force);
void cell_drift_gpart(struct cell *c, const struct engine *e, int force);
+void cell_drift_spart(struct cell *c, const struct engine *e, int force);
void cell_drift_multipole(struct cell *c, const struct engine *e);
void cell_drift_all_multipoles(struct cell *c, const struct engine *e);
void cell_check_timesteps(struct cell *c);
@@ -912,7 +922,8 @@ __attribute__((always_inline)) INLINE static int cell_can_split_pair_hydro_task(
/* Note that since tasks are create after a rebuild no need to take */
/* into account any part motion (i.e. dx_max == 0 here) */
return c->split &&
- (space_stretch * kernel_gamma * c->hydro.h_max < 0.5f * c->dmin);
+ (space_stretch * kernel_gamma * c->hydro.h_max < 0.5f * c->dmin) &&
+ (space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin);
}
/**
@@ -930,7 +941,8 @@ __attribute__((always_inline)) INLINE static int cell_can_split_self_hydro_task(
/* Note: No need for more checks here as all the sub-pairs and sub-self */
/* tasks will be created. So no need to check for h_max */
return c->split &&
- (space_stretch * kernel_gamma * c->hydro.h_max < 0.5f * c->dmin);
+ (space_stretch * kernel_gamma * c->hydro.h_max < 0.5f * c->dmin) &&
+ (space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin);
}
/**
diff --git a/src/runner.c b/src/runner.c
index 6585458064..c9ac839568 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1774,6 +1774,21 @@ void runner_do_drift_gpart(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_drift_gpart);
}
+/**
+ * @brief Drift all spart in a cell.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_drift_spart(struct runner *r, struct cell *c, int timer) {
+
+ TIMER_TIC;
+
+ cell_drift_spart(c, r->e, 0);
+
+ if (timer) TIMER_TOC(timer_drift_spart);
+}
/**
* @brief Perform the first half-kick on all the active particles in a cell.
*
@@ -3162,6 +3177,9 @@ void *runner_main(void *data) {
case task_type_drift_part:
runner_do_drift_part(r, ci, 1);
break;
+ case task_type_drift_spart:
+ runner_do_drift_spart(r, ci, 1);
+ break;
case task_type_drift_gpart:
runner_do_drift_gpart(r, ci, 1);
break;
diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h
index a8cd8a94b1..f1210acad8 100644
--- a/src/runner_doiact_stars.h
+++ b/src/runner_doiact_stars.h
@@ -1434,6 +1434,7 @@ void DOSUB_SELF1_STARS(struct runner *r, struct cell *ci, int gettimer) {
if (ci->nodeID != engine_rank)
error("This function should not be called on foreign cells");
#endif
+
/* Should we even bother? */
if (ci->hydro.count == 0 || ci->stars.count == 0 ||
!cell_is_active_stars(ci, r->e))
--
GitLab