From 17157f90083538ef71d1ac34723f64d983db37a5 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sat, 25 Feb 2017 14:08:02 +0000
Subject: [PATCH] Drift multipoles independently of the g-particles.
---
src/cell.c | 37 ++++++++++++++++++++++++++++++++++++-
src/cell.h | 6 +++++-
src/engine.c | 2 +-
src/multipole.h | 14 ++++++++++++++
src/runner.c | 12 +++++++++---
src/runner.h | 3 ++-
src/space.c | 3 +++
7 files changed, 70 insertions(+), 7 deletions(-)
diff --git a/src/cell.c b/src/cell.c
index c5f414817e..734b855851 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1326,7 +1326,7 @@ void cell_set_super(struct cell *c, struct cell *super) {
}
/**
- * @brief Recursively drifts all particles and g-particles in a cell hierarchy.
+ * @brief Recursively drifts particles of all kinds in a cell hierarchy.
*
* @param c The #cell.
* @param e The #engine (to get ti_current).
@@ -1430,6 +1430,41 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
c->ti_old = ti_current;
}
+/**
+ * @brief Recursively drifts all multipoles in a cell hierarchy.
+ *
+ * @param c The #cell.
+ * @param e The #engine (to get ti_current).
+ */
+void cell_drift_multipole(struct cell *c, const struct engine *e) {
+
+ const double timeBase = e->timeBase;
+ const integertime_t ti_old_multipole = c->ti_old_multipole;
+ const integertime_t ti_current = e->ti_current;
+
+ /* Drift from the last time the cell was drifted to the current time */
+ const double dt = (ti_current - ti_old_multipole) * timeBase;
+
+ /* Check that we are actually going to move forward. */
+ if (ti_current < ti_old_multipole) error("Attempt to drift to the past");
+
+ /* Are we not in a leaf ? */
+ if (c->split) {
+
+ /* Loop over the progeny and drift the multipoles. */
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) cell_drift_particles(c->progeny[k], e);
+
+ } else if (ti_current > ti_old_multipole) {
+
+ /* Drift the multipole */
+ multipole_drift(c->multipole, dt);
+ }
+
+ /* Update the time of the last drift */
+ c->ti_old_multipole = ti_current;
+}
+
/**
* @brief Recursively checks that all particles in a cell have a time-step
*/
diff --git a/src/cell.h b/src/cell.h
index 6c90f38c51..dc0e8886e0 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -230,9 +230,12 @@ struct cell {
/*! Maximum beginning of (integer) time step in this cell. */
integertime_t ti_beg_max;
- /*! Last (integer) time the cell's content was drifted forward in time. */
+ /*! Last (integer) time the cell's particle was drifted forward in time. */
integertime_t ti_old;
+ /*! Last (integer) time the cell's multipole was drifted forward in time. */
+ integertime_t ti_old_multipole;
+
/*! Minimum dimension, i.e. smallest edge of this cell (min(width)). */
float dmin;
@@ -343,6 +346,7 @@ int cell_is_drift_needed(struct cell *c, const struct engine *e);
int cell_unskip_tasks(struct cell *c, struct scheduler *s);
void cell_set_super(struct cell *c, struct cell *super);
void cell_drift_particles(struct cell *c, const struct engine *e);
+void cell_drift_multipole(struct cell *c, const struct engine *e);
void cell_check_timesteps(struct cell *c);
#endif /* SWIFT_CELL_H */
diff --git a/src/engine.c b/src/engine.c
index e46d0598f8..deae35f865 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3144,7 +3144,7 @@ void engine_unskip(struct engine *e) {
void engine_drift_all(struct engine *e) {
const ticks tic = getticks();
- threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells_top,
+ threadpool_map(&e->threadpool, runner_do_drift_all_mapper, e->s->cells_top,
e->s->nr_cells, sizeof(struct cell), 1, e);
#ifdef SWIFT_DEBUG_CHECKS
diff --git a/src/multipole.h b/src/multipole.h
index 69195f8ec7..9521bd51db 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -185,6 +185,20 @@ INLINE static int multipole_equal(const struct multipole *ma,
return 1;
}
+/**
+ * @brief Drifts a #multipole forward in time.
+ *
+ * @param m The #multipole.
+ * @param dt The drift time-step.
+ */
+INLINE static void multipole_drift(struct multipole *m, double dt) {
+
+ /* Move the whole thing according to bulk motion */
+ m->CoM[0] += m->vel[0];
+ m->CoM[1] += m->vel[1];
+ m->CoM[2] += m->vel[2];
+}
+
/**
* @brief Applies the forces due to particles j onto particles i directly.
*
diff --git a/src/runner.c b/src/runner.c
index 6f2f32c1f5..ff78e5d679 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -826,15 +826,21 @@ void runner_do_drift_particles(struct runner *r, struct cell *c, int timer) {
* @param num_elements Chunk size.
* @param extra_data Pointer to an #engine.
*/
-void runner_do_drift_mapper(void *map_data, int num_elements,
- void *extra_data) {
+void runner_do_drift_all_mapper(void *map_data, int num_elements,
+ void *extra_data) {
struct engine *e = (struct engine *)extra_data;
struct cell *cells = (struct cell *)map_data;
for (int ind = 0; ind < num_elements; ind++) {
struct cell *c = &cells[ind];
- if (c != NULL && c->nodeID == e->nodeID) cell_drift_particles(c, e);
+ if (c != NULL && c->nodeID == e->nodeID) {
+ /* Drift all the particles */
+ cell_drift_particles(c, e);
+
+ /* Drift the multipole */
+ if (e->policy & engine_policy_self_gravity) cell_drift_multipole(c, e);
+ }
}
}
diff --git a/src/runner.h b/src/runner.h
index 49ad5c9835..ec63eb3ec9 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -66,6 +66,7 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer);
void *runner_main(void *data);
void runner_do_unskip_mapper(void *map_data, int num_elements,
void *extra_data);
-void runner_do_drift_mapper(void *map_data, int num_elements, void *extra_data);
+void runner_do_drift_all_mapper(void *map_data, int num_elements,
+ void *extra_data);
#endif /* SWIFT_RUNNER_H */
diff --git a/src/space.c b/src/space.c
index 180b08c8f3..39a2ddbcd6 100644
--- a/src/space.c
+++ b/src/space.c
@@ -441,6 +441,7 @@ void space_regrid(struct space *s, int verbose) {
c->scount = 0;
c->super = c;
c->ti_old = ti_old;
+ c->ti_old_multipole = ti_old;
if (s->gravity) c->multipole = &s->multipoles_top[cid];
}
@@ -922,6 +923,7 @@ void space_rebuild(struct space *s, int verbose) {
for (int k = 0; k < s->nr_cells; k++) {
struct cell *restrict c = &cells_top[k];
c->ti_old = ti_old;
+ c->ti_old_multipole = ti_old;
c->parts = finger;
c->xparts = xfinger;
c->gparts = gfinger;
@@ -2030,6 +2032,7 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->gcount = 0;
cp->scount = 0;
cp->ti_old = c->ti_old;
+ cp->ti_old_multipole = c->ti_old_multipole;
cp->loc[0] = c->loc[0];
cp->loc[1] = c->loc[1];
cp->loc[2] = c->loc[2];
--
GitLab