diff --git a/src/cell.c b/src/cell.c index ccafae21f112d6504dbb4de8b6c915102d08858a..9cf64b89005c497d70feccb6bc4067bdbfb8ee02 100644 --- a/src/cell.c +++ b/src/cell.c @@ -262,7 +262,7 @@ int cell_pack_ti_ends(struct cell *c, int *ti_ends) { * @brief Unpack the time information of a given cell and its sub-cells. * * @param c The #cell - * @param s ti_ends The time information to unpack + * @param ti_ends The time information to unpack * * @return The number of cells created. */ @@ -794,9 +794,9 @@ void cell_check_multipole(struct cell *c, void *data) { } /** - * @brief Frees up the memory allocated for this #cell + * @brief Frees up the memory allocated for this #cell. * - * @param c The #cell + * @param c The #cell. */ void cell_clean(struct cell *c) { @@ -806,3 +806,34 @@ void cell_clean(struct cell *c) { for (int k = 0; k < 8; k++) if (c->progeny[k]) cell_clean(c->progeny[k]); } + +/** + * @brief Checks whether a given cell needs drifting or not. + * + * @param c the #cell. + * @param ti_current The current time on the integer time-line. + * + * @return 1 If the cell needs drifting, 0 otherwise. + */ +int cell_is_drift_needed(struct cell *c, int ti_current) { + + /* Do we have at least one active particle in the cell ?*/ + if (c->ti_end_min == ti_current) return 1; + + /* Loop over the pair tasks that involve this cell */ + for (struct link *l = c->density; l != NULL; l = l->next) { + + if (l->t->type != task_type_pair && l->t->type != task_type_sub_pair) + continue; + + /* Left or right? */ + if (l->t->ci == c) { + if (l->t->cj->ti_end_min == ti_current) return 1; + } else if (l->t->cj == c) { + if (l->t->ci->ti_end_min == ti_current) return 1; + } + } + + /* No neighbouring cell has active particles. Drift not necessary */ + return 0; +} diff --git a/src/cell.h b/src/cell.h index ee847e21cc4fe8a5859897ef2c2dc6dc9d367cb3..605d0c85444c1ac3b4147ef2fb23df52b2694ff5 100644 --- a/src/cell.h +++ b/src/cell.h @@ -47,6 +47,16 @@ struct space; /* Global variables. */ extern int cell_next_tag; +/* Mini struct to link cells to tasks. Used as a linked list. */ +struct link { + + /* The task pointer. */ + struct task *t; + + /* The next pointer. */ + struct link *next; +}; + /* Packed cell. */ struct pcell { @@ -79,6 +89,9 @@ struct cell { /* Minimum and maximum end of time step in this cell. */ int ti_end_min, ti_end_max; + /* Last time the cell's content was drifted forward in time. */ + int ti_old; + /* Minimum dimension, i.e. smallest edge of this cell. */ float dmin; @@ -208,5 +221,6 @@ int cell_are_neighbours(const struct cell *restrict ci, const struct cell *restrict cj); void cell_check_multipole(struct cell *c, void *data); void cell_clean(struct cell *c); +int cell_is_drift_needed(struct cell *c, int ti_current); #endif /* SWIFT_CELL_H */ diff --git a/src/engine.c b/src/engine.c index 42e7466f79a7e786cebd9cb5832f9f10e0c7641e..dd55fe25fd7938dc8a282c324de71fea80c8ca24 100644 --- a/src/engine.c +++ b/src/engine.c @@ -2088,6 +2088,14 @@ void engine_prepare(struct engine *e) { /* Did this not go through? */ if (rebuild) { + + /* First drift all particles to the current time */ + e->drift_all = 1; + threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells, + e->s->nr_cells, sizeof(struct cell), 1, e); + e->drift_all = 0; + + /* And now rebuild */ engine_rebuild(e); } @@ -2468,8 +2476,10 @@ void engine_step(struct engine *e) { snapshot_drift_time = e->timeStep; /* Drift everybody to the snapshot position */ + e->drift_all = 1; threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells, e->s->nr_cells, sizeof(struct cell), 1, e); + e->drift_all = 0; /* Dump... */ engine_dump_snapshot(e); @@ -2486,7 +2496,7 @@ void engine_step(struct engine *e) { e->timeOld = e->ti_old * e->timeBase + e->timeBegin; e->timeStep = (e->ti_current - e->ti_old) * e->timeBase + snapshot_drift_time; - /* Drift everybody */ + /* Drift only the necessary particles */ threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells, e->s->nr_cells, sizeof(struct cell), 1, e); @@ -2915,6 +2925,7 @@ void engine_init(struct engine *e, struct space *s, e->timeStep = 0.; e->timeBase = 0.; e->timeBase_inv = 0.; + e->drift_all = 0; e->internalUnits = internal_units; e->timeFirstSnapshot = parser_get_param_double(params, "Snapshots:time_first"); diff --git a/src/engine.h b/src/engine.h index e0c0f6a92e98a3c9a59f7db6a8ae930442ea5cac..61534ebfcde00e580f72a8d351280e157862728b 100644 --- a/src/engine.h +++ b/src/engine.h @@ -81,16 +81,6 @@ extern int engine_rank; /* The maximal number of timesteps in a simulation */ #define max_nr_timesteps (1 << 28) -/* Mini struct to link cells to density/force tasks. */ -struct link { - - /* The task pointer. */ - struct task *t; - - /* The next pointer. */ - struct link *next; -}; - /* Data structure for the engine. */ struct engine { @@ -139,6 +129,9 @@ struct engine { /* Minimal ti_end for the next time-step */ int ti_end_min; + /* Are we drifting all particles now ? */ + int drift_all; + /* Number of particles updated */ size_t updates, g_updates; diff --git a/src/runner.c b/src/runner.c index 23fd820d119230a499a87a447562fc631142655e..78491ec1dbf1173d8c9fd301dcb65199d086c10c 100644 --- a/src/runner.c +++ b/src/runner.c @@ -592,15 +592,19 @@ void runner_do_ghost(struct runner *r, struct cell *c) { static void runner_do_drift(struct cell *c, struct engine *e) { const double timeBase = e->timeBase; - const double dt = (e->ti_current - e->ti_old) * timeBase; - const int ti_old = e->ti_old; + const int ti_old = c->ti_old; const int ti_current = e->ti_current; - struct part *const parts = c->parts; struct xpart *const xparts = c->xparts; struct gpart *const gparts = c->gparts; - float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f; + /* Do we need to drift ? */ + if (!e->drift_all && !cell_is_drift_needed(c, ti_current)) return; + + /* Drift from the last time the cell was drifted to the current time */ + const double dt = (ti_current - ti_old) * timeBase; + + float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f; double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0; double mom[3] = {0.0, 0.0, 0.0}; double ang_mom[3] = {0.0, 0.0, 0.0}; @@ -722,6 +726,9 @@ static void runner_do_drift(struct cell *c, struct engine *e) { c->ang_mom[0] = ang_mom[0]; c->ang_mom[1] = ang_mom[1]; c->ang_mom[2] = ang_mom[2]; + + /* Update the time of the last drift */ + c->ti_old = ti_current; } /** diff --git a/src/space.c b/src/space.c index 57b7279a508f10e3647105562e5974072a0c57f0..def8adad3d023efe5daf552328d523d3cfc37429 100644 --- a/src/space.c +++ b/src/space.c @@ -164,6 +164,7 @@ void space_regrid(struct space *s, double cell_max, int verbose) { const size_t nr_parts = s->nr_parts; struct cell *restrict c; ticks tic = getticks(); + const int ti_current = (s->e != NULL) ? s->e->ti_current : 0; /* Run through the cells and get the current h_max. */ // tic = getticks(); @@ -294,6 +295,7 @@ void space_regrid(struct space *s, double cell_max, int verbose) { c->count = 0; c->gcount = 0; c->super = c; + c->ti_old = ti_current; lock_init(&c->lock); } @@ -392,6 +394,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) { size_t nr_parts = s->nr_parts; size_t nr_gparts = s->nr_gparts; struct cell *restrict cells = s->cells; + const int ti_current = (s->e != NULL) ? s->e->ti_current : 0; const double ih[3] = {s->iwidth[0], s->iwidth[1], s->iwidth[2]}; const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]}; @@ -643,6 +646,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) { struct gpart *gfinger = s->gparts; for (int k = 0; k < s->nr_cells; k++) { struct cell *restrict c = &cells[k]; + c->ti_old = ti_current; c->parts = finger; c->xparts = xfinger; c->gparts = gfinger; @@ -1234,6 +1238,7 @@ void space_split_mapper(void *map_data, int num_elements, void *extra_data) { /* Unpack the inputs. */ struct space *s = (struct space *)extra_data; struct cell *cells = (struct cell *)map_data; + struct engine *e = s->e; for (int ind = 0; ind < num_elements; ind++) { @@ -1263,6 +1268,7 @@ void space_split_mapper(void *map_data, int num_elements, void *extra_data) { temp = space_getcell(s); temp->count = 0; temp->gcount = 0; + temp->ti_old = e->ti_current; temp->loc[0] = c->loc[0]; temp->loc[1] = c->loc[1]; temp->loc[2] = c->loc[2];