Commit ffb2e88d authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Drift before init for active cells

parent 0b6ee946
......@@ -721,7 +721,8 @@ void cell_check_drift_point(struct cell *c, void *data) {
/* for (int i = 0; i < c->count; ++i) */
/* if (c->parts[i].ti_old != ti_current) */
/* error( */
/* "Particle in an incorrect time-zone! part->ti_old=%d c->ti_old=%d " */
/* "Particle in an incorrect time-zone! part->ti_old=%d c->ti_old=%d "
*/
/* "ti_current=%d", */
/* c->parts[i].ti_old, c->ti_old, ti_current); */
}
......@@ -1009,11 +1010,11 @@ void cell_drift(struct cell *c, const struct engine *e) {
const double dt = (ti_current - ti_old) * timeBase;
float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
//message("DRFIT ! ti_old=%d ti_current=%d", ti_old, ti_current);
// message("DRFIT ! ti_old=%d ti_current=%d", ti_old, ti_current);
/* Check that we are actually going to move forward. */
if (ti_current < ti_old) error("Attempt to drift to the past");
//if (ti_current == ti_old) return;
// if (ti_current == ti_old) return;
/* Are we not in a leaf ? */
if (c->split) {
......@@ -1057,7 +1058,7 @@ void cell_drift(struct cell *c, const struct engine *e) {
/* Drift... */
drift_part(p, xp, dt, timeBase, ti_old, ti_current);
//p->ti_old = ti_current;
// p->ti_old = ti_current;
/* Compute (square of) motion since last cell construction */
const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
......
......@@ -142,6 +142,11 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0, 0, c,
NULL, 0);
c->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0, 0,
c, NULL, 0);
scheduler_addunlock(s, c->drift, c->init);
/* Generate the ghost task. */
if (is_hydro)
c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0,
......@@ -1987,10 +1992,10 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
else if (t->type == task_type_self || t->type == task_type_sub_self) {
/* Local pointers. */
const struct cell *ci = t->ci;
// const struct cell *ci = t->ci;
/* Activate the drift */
if (ci->drift) scheduler_activate(s, ci->drift);
// if (ci->drift) scheduler_activate(s, ci->drift);
}
/* Pair? */
......@@ -2001,8 +2006,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
const struct cell *cj = t->cj;
/* Activate the drift on both sides */
if (ci->drift) scheduler_activate(s, ci->drift);
if (cj->drift) scheduler_activate(s, cj->drift);
// if (ci->drift) scheduler_activate(s, ci->drift);
// if (cj->drift) scheduler_activate(s, cj->drift);
/* Too much particle movement? */
if (t->tight &&
......@@ -2100,6 +2105,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (t->ci->ti_end_min <= ti_end) scheduler_activate(s, t);
}
/* Drift? */
else if (t->type == task_type_drift) {
if (t->ci->ti_end_min <= ti_end) scheduler_activate(s, t);
}
/* Init? */
else if (t->type == task_type_init) {
if (t->ci->ti_end_min <= ti_end) scheduler_activate(s, t);
......@@ -2482,7 +2492,8 @@ void engine_skip_force_and_kick(struct engine *e) {
/* Skip everything that updates the particles */
if (t->subtype == task_subtype_force || t->type == task_type_kick ||
t->type == task_type_cooling || t->type == task_type_sourceterms)
t->type == task_type_cooling || t->type == task_type_sourceterms ||
t->type == task_type_drift)
t->skip = 1;
}
}
......
......@@ -212,7 +212,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part *restrict p, struct xpart *restrict xp) {
//p->ti_old = 0;
// p->ti_old = 0;
p->ti_begin = 0;
p->ti_end = 0;
xp->v_full[0] = p->v[0];
......
......@@ -55,7 +55,7 @@ struct part {
/* Pointer to corresponding gravity part. */
struct gpart* gpart;
/* Particle position. */
double x[3];
......@@ -77,8 +77,8 @@ struct part {
/* Particle time of end of time-step. */
int ti_end;
//int ti_old;
// int ti_old;
/* Particle density. */
float rho;
......@@ -88,7 +88,7 @@ struct part {
/* Entropy time derivative */
float entropy_dt;
//union {
union {
struct {
......@@ -130,7 +130,7 @@ struct part {
float h_dt;
} force;
//};
};
} SWIFT_STRUCT_ALIGN;
......
......@@ -736,11 +736,11 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
if (count) {
message("Smoothing length failed to converge on %i particles.", count);
for(int i=0; i<count; ++i) {
struct part *restrict p = &parts[pid[i]];
for (int i = 0; i < count; ++i) {
struct part *restrict p = &parts[pid[i]];
struct xpart *restrict xp = &xparts[pid[i]];
printParticle_single(p, xp);
printParticle_single(p, xp);
}
}
......
......@@ -721,8 +721,8 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
/* Anything to do here? */
if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
if(!cell_is_drifted(ci, e)) cell_drift(ci, e);
if(!cell_is_drifted(cj, e)) cell_drift(cj, e);
if (!cell_is_drifted(ci, e)) cell_drift(ci, e);
if (!cell_is_drifted(cj, e)) cell_drift(cj, e);
/* Get the sort ID. */
double shift[3] = {0.0, 0.0, 0.0};
......@@ -762,8 +762,6 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
const double di = sort_i[pid].d + hi * kernel_gamma + dx_max - rshift;
if (di < dj_min) continue;
double pix[3];
for (int k = 0; k < 3; k++) pix[k] = pi->x[k] - shift[k];
const float hig2 = hi * hi * kernel_gamma2;
......@@ -774,7 +772,6 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
/* Get a pointer to the jth particle. */
struct part *restrict pj = &parts_j[sort_j[pjd].i];
/* Compute the pairwise distance. */
float r2 = 0.0f;
float dx[3];
......@@ -827,8 +824,6 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
const double dj = sort_j[pjd].d - hj * kernel_gamma - dx_max - rshift;
if (dj > di_max) continue;
double pjx[3];
for (int k = 0; k < 3; k++) pjx[k] = pj->x[k] + shift[k];
const float hjg2 = hj * hj * kernel_gamma2;
......@@ -839,8 +834,6 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
/* Get a pointer to the jth particle. */
struct part *restrict pi = &parts_i[sort_i[pid].i];
/* Compute the pairwise distance. */
float r2 = 0.0f;
float dx[3];
......@@ -923,9 +916,9 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
/* Anything to do here? */
if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
if(!cell_is_drifted(ci, e)) error("Cell ci not drifted");
if(!cell_is_drifted(cj, e)) error("Cell cj not drifted");
if (!cell_is_drifted(ci, e)) error("Cell ci not drifted");
if (!cell_is_drifted(cj, e)) error("Cell cj not drifted");
/* Get the shift ID. */
double shift[3] = {0.0, 0.0, 0.0};
const int sid = space_getsid(e->s, &ci, &cj, shift);
......@@ -993,8 +986,6 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
const double di = sort_i[pid].d + hi * kernel_gamma + dx_max - rshift;
if (di < dj_min) continue;
double pix[3];
for (int k = 0; k < 3; k++) pix[k] = pi->x[k] - shift[k];
const float hig2 = hi * hi * kernel_gamma2;
......@@ -1009,8 +1000,6 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
struct part *restrict pj = &parts_j[sortdt_j[pjd].i];
const float hj = pj->h;
/* Compute the pairwise distance. */
float r2 = 0.0f;
float dx[3];
......@@ -1062,8 +1051,6 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
struct part *restrict pj = &parts_j[sort_j[pjd].i];
const float hj = pj->h;
/* Compute the pairwise distance. */
float r2 = 0.0f;
float dx[3];
......@@ -1143,8 +1130,6 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
const double dj = sort_j[pjd].d - hj * kernel_gamma - dx_max - rshift;
if (dj > di_max) continue;
double pjx[3];
for (int k = 0; k < 3; k++) pjx[k] = pj->x[k] + shift[k];
const float hjg2 = hj * hj * kernel_gamma2;
......@@ -1159,8 +1144,6 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
struct part *restrict pi = &parts_i[sortdt_i[pid].i];
const float hi = pi->h;
/* Compute the pairwise distance. */
float r2 = 0.0f;
float dx[3];
......@@ -1211,8 +1194,6 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
struct part *restrict pi = &parts_i[sort_i[pid].i];
const float hi = pi->h;
/* Compute the pairwise distance. */
float r2 = 0.0f;
float dx[3];
......@@ -1328,7 +1309,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
if (!cell_is_active(c, e)) return;
if(!cell_is_drifted(c, e)) cell_drift(c, e);
if (!cell_is_drifted(c, e)) cell_drift(c, e);
struct part *restrict parts = c->parts;
const int count = c->count;
......@@ -1357,8 +1338,6 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
const float hi = pi->h;
const float hig2 = hi * hi * kernel_gamma2;
/* Is the ith particle inactive? */
if (!part_is_active(pi, e)) {
......@@ -1369,8 +1348,6 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
struct part *restrict pj = &parts[indt[pjd]];
const float hj = pj->h;
/* Compute the pairwise distance. */
float r2 = 0.0f;
float dx[3];
......@@ -1425,8 +1402,6 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
struct part *restrict pj = &parts[pjd];
const float hj = pj->h;
/* Compute the pairwise distance. */
float r2 = 0.0f;
float dx[3];
......@@ -1567,7 +1542,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
if (!cell_is_active(c, e)) return;
if(!cell_is_drifted(c, e)) error("Cell is not drifted");
if (!cell_is_drifted(c, e)) error("Cell is not drifted");
struct part *restrict parts = c->parts;
const int count = c->count;
......@@ -1596,8 +1571,6 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
const float hi = pi->h;
const float hig2 = hi * hi * kernel_gamma2;
/* Is the ith particle not active? */
if (!part_is_active(pi, e)) {
......@@ -1608,8 +1581,6 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
struct part *restrict pj = &parts[indt[pjd]];
const float hj = pj->h;
/* Compute the pairwise distance. */
float r2 = 0.0f;
float dx[3];
......@@ -1664,8 +1635,6 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
struct part *restrict pj = &parts[pjd];
const float hj = pj->h;
/* Compute the pairwise distance. */
float r2 = 0.0f;
float dx[3];
......
......@@ -216,7 +216,6 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
to make sure we get ci and cj swapped if needed. */
double shift[3];
const int sid = space_getsid(s->space, &ci, &cj, shift);
const int did = space_getdid(s->space, ci, cj);
/* Should this task be split-up? */
if (ci->split && cj->split &&
......@@ -617,28 +616,10 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
1 << sid, 0, ci, NULL, 0);
else
ci->sorts->flags |= (1 << sid);
lock_unlock_blind(&ci->lock);
scheduler_addunlock(s, ci->sorts, t);
/* Create the drift for ci. */
if (ci->drift == NULL) {
// ci->drift = scheduler_addtask(s, task_type_drift,
// task_subtype_none,
// 1 << sid, 0, ci, NULL, 0);
// scheduler_addunlock(s, ci->drift, ci->sorts);
// scheduler_addunlock(s, ci->drift, t);
}
/* if(!(ci->drift->flags & (1 << sid))) { */
/* scheduler_addunlock(s, ci->drift, ci->sorts); */
/* ci->drift->flags |= (1 << sid); */
/* } */
if (did == 0 || did > 31) message("did=%d 1<<did=%d", did, 1 << did);
lock_unlock_blind(&ci->lock);
/* Create the sort for cj. */
lock_lock(&cj->lock);
if (cj->sorts == NULL)
......@@ -646,26 +627,9 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
1 << sid, 0, cj, NULL, 0);
else
cj->sorts->flags |= (1 << sid);
lock_unlock_blind(&cj->lock);
scheduler_addunlock(s, cj->sorts, t);
/* Create the drift for cj. */
if (cj->drift == NULL) {
// cj->drift = scheduler_addtask(s, task_type_drift,
// task_subtype_none,
// 1 << sid, 0, cj, NULL, 0);
// scheduler_addunlock(s, cj->drift, cj->sorts);
// scheduler_addunlock(s, cj->drift, t);
// scheduler_addunlock(s, cj->drift, cj->init);
}
/* if(!(cj->drift->flags & (1 << sid))) { */
/* scheduler_addunlock(s, cj->drift, cj->sorts); */
/* cj->drift->flags |= (1 << sid); */
/* } */
lock_unlock_blind(&cj->lock);
}
} /* pair interaction? */
......
......@@ -168,32 +168,6 @@ int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
return sid;
}
int space_getdid(struct space *s, struct cell *ci, struct cell *cj) {
/* Get the relative distance between the pairs, wrapping. */
const int periodic = s->periodic;
double dx[3];
double shift[3];
for (int k = 0; k < 3; k++) {
dx[k] = cj->loc[k] - ci->loc[k];
if (periodic && dx[k] < -s->dim[k] / 2)
shift[k] = s->dim[k];
else if (periodic && dx[k] > s->dim[k] / 2)
shift[k] = -s->dim[k];
else
shift[k] = 0.0;
dx[k] += shift[k];
}
/* Get the drift index. */
int did = 0;
for (int k = 0; k < 3; k++)
did = 3 * did + ((dx[k] < 0.0) ? 0 : ((dx[k] > 0.0) ? 2 : 1));
/* Return the drift ID. */
return did;
}
/**
* @brief Recursively dismantle a cell tree.
*
......
......@@ -150,7 +150,6 @@ void space_gparts_sort(struct space *s, int *ind, size_t N, int min, int max,
struct cell *space_getcell(struct space *s);
int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
double *shift);
int space_getdid(struct space *s, struct cell *ci, struct cell *cj);
void space_init(struct space *s, const struct swift_params *params,
double dim[3], struct part *parts, struct gpart *gparts,
size_t Npart, size_t Ngpart, int periodic, int gravity,
......
Supports Markdown
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