Commit 0058b193 authored by Matthieu Schaller's avatar Matthieu Schaller

Increase accuracy of drift. Synchronize part<->gpart positions when drifiting all.

parent 65e95435
......@@ -39,7 +39,7 @@
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_gpart(
struct gpart *restrict gp, float dt, double timeBase, integertime_t ti_old,
struct gpart *restrict gp, double dt, double timeBase, integertime_t ti_old,
integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
......@@ -75,7 +75,7 @@ __attribute__((always_inline)) INLINE static void drift_gpart(
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_part(
struct part *restrict p, struct xpart *restrict xp, float dt,
struct part *restrict p, struct xpart *restrict xp, double dt,
double timeBase, integertime_t ti_old, integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
......@@ -119,7 +119,7 @@ __attribute__((always_inline)) INLINE static void drift_part(
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_spart(
struct spart *restrict sp, float dt, double timeBase, integertime_t ti_old,
struct spart *restrict sp, double dt, double timeBase, integertime_t ti_old,
integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
......
......@@ -1872,6 +1872,14 @@ void engine_count_and_link_tasks(struct engine *e) {
if (finger != NULL) scheduler_addunlock(sched, t, finger->drift_part);
}
/* Link drift tasks to the next-higher drift task. */
else if (t->type == task_type_drift_gpart) {
struct cell *finger = ci->parent;
while (finger != NULL && finger->drift_gpart == NULL)
finger = finger->parent;
if (finger != NULL) scheduler_addunlock(sched, t, finger->drift_gpart);
}
/* Link self tasks to cells. */
else if (t->type == task_type_self) {
atomic_inc(&ci->nr_tasks);
......@@ -3055,7 +3063,7 @@ void engine_print_stats(struct engine *e) {
e->policy & engine_policy_self_gravity);
/* Be verbose about this */
message("Saving statistics at t=%e.", e->time);
if(e->nodeID == 0) message("Saving statistics at t=%e.", e->time);
#else
if (e->verbose) message("Saving statistics at t=%e.", e->time);
#endif
......@@ -3535,10 +3543,15 @@ void engine_drift_all(struct engine *e) {
threadpool_map(&e->threadpool, engine_do_drift_all_mapper, e->s->cells_top,
e->s->nr_cells, sizeof(struct cell), 1, e);
/* Synchronize particle positions */
space_synchronize_particle_positions(e->s);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all cells have been drifted to the current time. */
space_check_drift_point(e->s, e->ti_current,
e->policy & engine_policy_self_gravity);
part_verify_links(e->s->parts, e->s->gparts, e->s->sparts, e->s->nr_parts,
e->s->nr_gparts, e->s->nr_sparts, e->verbose);
#endif
if (e->verbose)
......@@ -3837,7 +3850,7 @@ void engine_dump_snapshot(struct engine *e) {
e->policy & engine_policy_self_gravity);
/* Be verbose about this */
message("writing snapshot at t=%e.", e->time);
if (e->nodeID == 0) message("writing snapshot at t=%e.", e->time);
#else
if (e->verbose) message("writing snapshot at t=%e.", e->time);
#endif
......
......@@ -2431,6 +2431,58 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
}
}
void space_synchronize_particle_positions_mapper(void *map_data, int nr_gparts,
void *extra_data) {
/* Unpack the data */
struct gpart *restrict gparts = (struct gpart *)map_data;
struct space *s = (struct space*) extra_data;
for (int k = 0; k < nr_gparts; k++) {
/* Get the particle */
const struct gpart *restrict gp = &gparts[k];
if(gp->type == swift_type_dark_matter) continue;
else if(gp->type == swift_type_gas) {
/* Get it's gassy friend */
struct part *p = &s->parts[-gp->id_or_neg_offset];
struct xpart *xp = &s->xparts[-gp->id_or_neg_offset];
/* Synchronize positions and velocities */
p->x[0] = gp->x[0];
p->x[1] = gp->x[1];
p->x[2] = gp->x[2];
xp->v_full[0] = gp->v_full[0];
xp->v_full[1] = gp->v_full[1];
xp->v_full[2] = gp->v_full[2];
}
else if(gp->type == swift_type_star) {
/* Get it's stellar friend */
struct spart *sp = &s->sparts[-gp->id_or_neg_offset];
/* Synchronize positions */
sp->x[0] = gp->x[0];
sp->x[1] = gp->x[1];
sp->x[2] = gp->x[2];
}
}
}
void space_synchronize_particle_positions(struct space *s) {
if ((s->nr_gparts > 0 && s->nr_parts > 0) ||
(s->nr_gparts > 0 && s->nr_sparts > 0))
threadpool_map(&s->e->threadpool, space_synchronize_particle_positions_mapper,
s->gparts, s->nr_gparts, sizeof(struct gpart), 1000, (void*)s);
}
/**
* @brief Initialises all the particles by setting them into a valid state
*
......
......@@ -209,6 +209,7 @@ void space_gparts_get_cell_index(struct space *s, int *gind, struct cell *cells,
int verbose);
void space_sparts_get_cell_index(struct space *s, int *sind, struct cell *cells,
int verbose);
void space_synchronize_particle_positions(struct space *s);
void space_do_parts_sort();
void space_do_gparts_sort();
void space_do_sparts_sort();
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment