diff --git a/src/cell.c b/src/cell.c index 9bd2a104c90e3d6f7c7a7b71932e741677cd6f16..5557dd15d9cb9bf806e0f7ee146cf098a711a850 100644 --- a/src/cell.c +++ b/src/cell.c @@ -49,6 +49,7 @@ /* Local headers. */ #include "active.h" #include "atomic.h" +#include "drift.h" #include "error.h" #include "gravity.h" #include "hydro.h" @@ -985,3 +986,109 @@ void cell_set_super(struct cell *c, struct cell *super) { for (int k = 0; k < 8; k++) if (c->progeny[k] != NULL) cell_set_super(c->progeny[k], super); } + +void cell_drift(struct cell *c, struct engine *e) { + + const double timeBase = e->timeBase; + 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; + + /* 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; + + /* Check that we are actually going to move forward. */ + if (ti_current <= ti_old) return; + + /* Are we in a leaf ? */ + if (c->split) { + + /* 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]; + cell_drift(cp, e); + dx_max = max(dx_max, cp->dx_max); + h_max = max(h_max, cp->h_max); + } + + } else { + + /* Loop over all the g-particles in the cell */ + const size_t nr_gparts = c->gcount; + for (size_t k = 0; k < nr_gparts; k++) { + + /* Get a handle on the gpart. */ + struct gpart *const gp = &gparts[k]; + + /* Drift... */ + drift_gpart(gp, dt, timeBase, ti_old, ti_current); + + /* Compute (square of) motion since last cell construction */ + const float dx2 = gp->x_diff[0] * gp->x_diff[0] + + gp->x_diff[1] * gp->x_diff[1] + + gp->x_diff[2] * gp->x_diff[2]; + dx2_max = (dx2_max > dx2) ? dx2_max : dx2; + } + + /* Loop over all the particles in the cell */ + const size_t nr_parts = c->count; + for (size_t k = 0; k < nr_parts; k++) { + + /* Get a handle on the part. */ + struct part *const p = &parts[k]; + struct xpart *const xp = &xparts[k]; + + /* Drift... */ + drift_part(p, xp, dt, timeBase, ti_old, ti_current); + + /* Compute (square of) motion since last cell construction */ + const float dx2 = xp->x_diff[0] * xp->x_diff[0] + + xp->x_diff[1] * xp->x_diff[1] + + xp->x_diff[2] * xp->x_diff[2]; + dx2_max = (dx2_max > dx2) ? dx2_max : dx2; + + /* Maximal smoothing length */ + h_max = (h_max > p->h) ? h_max : p->h; + } + + /* Now, get the maximal particle motion from its square */ + dx_max = sqrtf(dx2_max); + + /* Set ti_old on the sub-cells */ + cell_set_ti_old(c, e->ti_current); + + } /* Check that we are actually going to move forward. */ + + /* Store the values */ + c->h_max = h_max; + c->dx_max = dx_max; + + /* Update the time of the last drift */ + c->ti_old = ti_current; +} + +/** + * Set ti_old of a #cell and all its progenies to a new value. + * + * @param c The #cell. + * @param ti_current The new value of ti_old. + */ +void cell_set_ti_old(struct cell *c, int ti_current) { + + /* Set this cell */ + c->ti_old = ti_current; + + /* Recurse */ + if (c->split) { + for (int k = 0; k < 8; ++k) { + if (c->progeny[k] != NULL) { + struct cell *cp = c->progeny[k]; + cell_set_ti_old(cp, ti_current); + } + } + } +} diff --git a/src/cell.h b/src/cell.h index d978de71f7fae8eeb73b527b064b6d88db6ccc1a..8a7301a57ad60a6ae8cbae755ced15aa8beeb06f 100644 --- a/src/cell.h +++ b/src/cell.h @@ -298,5 +298,6 @@ void cell_check_drift_point(struct cell *c, void *data); 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(struct cell *c, struct engine *e); +void cell_set_ti_old(struct cell *c, int ti_current); #endif /* SWIFT_CELL_H */ diff --git a/src/engine.c b/src/engine.c index be58038f5f544fc6a12dbf5137be6ff323f627db..d6a1dcf8b3d80cd816ec5971722bb57555ddd900 100644 --- a/src/engine.c +++ b/src/engine.c @@ -2702,7 +2702,7 @@ void engine_drift_all(struct engine *e) { e->drift_all = 1; const ticks tic = getticks(); - threadpool_map(&e->threadpool, runner_do_unskip_mapper, e->s->cells_top, + threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells_top, e->s->nr_cells, sizeof(struct cell), 1, e); e->drift_all = e->policy & engine_policy_drift_all; diff --git a/src/runner.c b/src/runner.c index fa42af4e99d44351eeb4f70a6545056c23ff3d93..334aa556048739202fdbc34458d917be3a9e9b93 100644 --- a/src/runner.c +++ b/src/runner.c @@ -779,23 +779,23 @@ static void runner_do_unskip(struct cell *c, struct engine *e, int drift) { /* Now, we can drift */ /* Get some information first */ - const double timeBase = e->timeBase; + // const double timeBase = e->timeBase; 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; + // struct part *const parts = c->parts; + // struct xpart *const xparts = c->xparts; + // struct gpart *const gparts = c->gparts; /* 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; + // const double dt = (ti_current - ti_old) * timeBase; + // float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f; /* No children? */ if (!c->split) { /* Check that we are actually going to move forward. */ if (ti_current > ti_old) { - +#if 0 /* Loop over all the g-particles in the cell */ const size_t nr_gparts = c->gcount; for (size_t k = 0; k < nr_gparts; k++) { @@ -836,13 +836,15 @@ static void runner_do_unskip(struct cell *c, struct engine *e, int drift) { /* Now, get the maximal particle motion from its square */ dx_max = sqrtf(dx2_max); - +#endif } /* Check that we are actually going to move forward. */ else { +#if 0 /* ti_old == ti_current, just keep the current cell values. */ h_max = c->h_max; dx_max = c->dx_max; +#endif } } @@ -853,20 +855,23 @@ static void runner_do_unskip(struct cell *c, struct engine *e, int drift) { for (int k = 0; k < 8; k++) if (c->progeny[k] != NULL) { struct cell *cp = c->progeny[k]; - + message("aaa"); /* Recurse. */ runner_do_unskip(cp, e, drift); +#if 0 dx_max = max(dx_max, cp->dx_max); h_max = max(h_max, cp->h_max); +#endif } } - +#if 0 /* Store the values */ c->h_max = h_max; c->dx_max = dx_max; /* Update the time of the last drift */ c->ti_old = ti_current; +#endif } /** @@ -876,7 +881,6 @@ static void runner_do_unskip(struct cell *c, struct engine *e, int drift) { * @param num_elements Chunk size. * @param extra_data Pointer to an #engine. */ - void runner_do_unskip_mapper(void *map_data, int num_elements, void *extra_data) { @@ -893,7 +897,22 @@ void runner_do_unskip_mapper(void *map_data, int num_elements, } } -void runner_do_drift(struct runner *r, struct cell *c, int timer) {} +void runner_do_drift(struct runner *r, struct cell *c, int timer) { + + cell_drift(c, r->e); +} + +void runner_do_drift_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(c, e); + } +} /** * @brief Kick particles in momentum space and collect statistics (floating diff --git a/src/runner.h b/src/runner.h index 5c89c97c6a00575ecd5966f007f92f27cc55f7a8..6aab8469a1eb9aba22dac31553424bbcce8f183f 100644 --- a/src/runner.h +++ b/src/runner.h @@ -58,5 +58,6 @@ 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); #endif /* SWIFT_RUNNER_H */