Skip to content

GitLab

Commit 0b6ee946 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

drift on demand working

parent 2d03f5ad
Hide whitespace changes
Inline Side-by-side
examples/SedovBlast_3D/sedov.yml
View file @ 0b6ee946
......@@ -6,6 +6,9 @@ InternalUnitSystem:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
#Scheduler:
# max_top_level_cells: 3
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
......
src/active.h
View file @ 0b6ee946
......@@ -36,7 +36,7 @@
* @param c The #cell.
* @param e The #engine containing information about the current time.
*/
__attribute__((always_inline)) INLINE static void cell_is_drifted(
__attribute__((always_inline)) INLINE static int cell_is_drifted(
const struct cell *c, const struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
......@@ -45,14 +45,35 @@ __attribute__((always_inline)) INLINE static void cell_is_drifted(
"Cell has been drifted too far forward in time! c->ti_old=%d "
"e->ti_current=%d",
c->ti_old, e->ti_current);
if (c->ti_old != e->ti_current) {
error(
"Cell has not been drifted to the current time c->ti_old=%d, "
"e->ti_current=%d",
c->ti_old, e->ti_current);
}
#endif
/* if (c->ti_old != e->ti_current) { */
/* int wrong = 0; */
/* for (int i = 0; i < c->count; ++i) { */
/* if (c->parts[i].ti_old < e->ti_current) ++wrong; */
/* } */
/* message( */
/* "Cell has not been drifted to the current time c->ti_old=%d, " */
/* "e->ti_current=%d wrong=%d c->count=%d c->drift=%p, c->depth=%d,
* c=%p, c->super=%p, c->parent=%p ", */
/* c->ti_old, e->ti_current, wrong, c->count, c->drift, c->depth, c,
* c->super, c->parent); */
/* cell_drift((struct cell*)c, e); */
/* message( */
/* "Cell has not been drifted to the current time c->ti_old=%d, " */
/* "e->ti_current=%d wrong=%d c->count=%d c->drift=%p, c->depth=%d,
* c=%p, c->super=%p, c->parent=%p ", */
/* c->ti_old, e->ti_current, wrong, c->count, c->drift, c->depth, c,
* c->super, c->parent); */
/* error("AAAAA"); */
/* } */
return (c->ti_old == e->ti_current);
}
/**
......
src/cell.c
View file @ 0b6ee946
......@@ -717,6 +717,13 @@ void cell_check_drift_point(struct cell *c, void *data) {
if (c->ti_old != ti_current)
error("Cell in an incorrect time-zone! c->ti_old=%d ti_current=%d",
c->ti_old, ti_current);
/* 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 " */
/* "ti_current=%d", */
/* c->parts[i].ti_old, c->ti_old, ti_current); */
}
/**
......@@ -880,8 +887,10 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
}
/* Activate the drift on both sides */
if (ci == c && cj != NULL && cj->drift != NULL) scheduler_activate(s, cj->drift);
if (cj == c && ci != NULL && ci->drift != NULL) scheduler_activate(s, ci->drift);
if (ci == c && cj != NULL && cj->drift != NULL)
scheduler_activate(s, cj->drift);
if (cj == c && ci != NULL && ci->drift != NULL)
scheduler_activate(s, ci->drift);
/* Check whether there was too much particle motion */
if (t->type == task_type_pair || t->type == task_type_sub_pair) {
......@@ -1000,8 +1009,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);
/* Check that we are actually going to move forward. */
if (ti_current <= ti_old) return;
if (ti_current < ti_old) error("Attempt to drift to the past");
//if (ti_current == ti_old) return;
/* Are we not in a leaf ? */
if (c->split) {
......@@ -1015,7 +1027,7 @@ void cell_drift(struct cell *c, const struct engine *e) {
h_max = max(h_max, cp->h_max);
}
} else {
} else if (ti_current >= ti_old) {
/* Loop over all the g-particles in the cell */
const size_t nr_gparts = c->gcount;
......@@ -1045,7 +1057,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] +
......@@ -1061,9 +1073,14 @@ void cell_drift(struct cell *c, const struct engine *e) {
dx_max = sqrtf(dx2_max);
/* Set ti_old on the sub-cells */
//cell_set_ti_old(c, e->ti_current);
cell_set_ti_old(c, e->ti_current);
} /* Check that we are actually going to move forward. */
else {
h_max = c->h_max;
dx_max = c->dx_max;
}
/* Store the values */
c->h_max = h_max;
......
src/engine.c
View file @ 0b6ee946
......@@ -2262,6 +2262,11 @@ void engine_prepare(struct engine *e, int nodrift) {
#endif
engine_rebuild(e);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all cells have been drifted to the current time */
space_check_drift_point(e->s, e->ti_current);
#endif
}
/* Re-rank the tasks every now and then. */
......@@ -2605,7 +2610,7 @@ void engine_step(struct engine *e) {
snapshot_drift_time = e->timeStep;
/* Drift everybody to the snapshot position */
e->drift_all = 1;
// e->drift_all = 1;
engine_drift_all(e);
/* Restore the default drifting policy */
......@@ -2629,8 +2634,9 @@ void engine_step(struct engine *e) {
if (e->nodeID == 0) {
/* Print some information to the screen */
printf(" %6d %14e %14e %10zu %10zu %21.3f\n", e->step, e->time,
e->timeStep, e->updates, e->g_updates, e->wallclock_time);
printf(" %6d %14e %d %14e %10zu %10zu %21.3f\n", e->step, e->time,
e->ti_current, e->timeStep, e->updates, e->g_updates,
e->wallclock_time);
fflush(stdout);
fprintf(e->file_timesteps, " %6d %14e %14e %10zu %10zu %21.3f\n", e->step,
......
src/hydro/Gadget2/hydro.h
View file @ 0b6ee946
......@@ -212,6 +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_begin = 0;
p->ti_end = 0;
xp->v_full[0] = p->v[0];
......
src/hydro/Gadget2/hydro_part.h
View file @ 0b6ee946
......@@ -50,6 +50,12 @@ struct xpart {
/* Data of a single particle. */
struct part {
/* Particle ID. */
long long id;
/* Pointer to corresponding gravity part. */
struct gpart* gpart;
/* Particle position. */
double x[3];
......@@ -71,6 +77,8 @@ struct part {
/* Particle time of end of time-step. */
int ti_end;
//int ti_old;
/* Particle density. */
float rho;
......@@ -80,7 +88,7 @@ struct part {
/* Entropy time derivative */
float entropy_dt;
union {
//union {
struct {
......@@ -122,13 +130,7 @@ struct part {
float h_dt;
} force;
};
/* Particle ID. */
long long id;
/* Pointer to corresponding gravity part. */
struct gpart* gpart;
//};
} SWIFT_STRUCT_ALIGN;
......
src/runner.c
View file @ 0b6ee946
......@@ -733,9 +733,17 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
}
}
if (count)
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]];
struct xpart *restrict xp = &xparts[pid[i]];
printParticle_single(p, xp);
}
}
/* Be clean */
free(pid);
}
......@@ -855,9 +863,9 @@ 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");
// message("aaa");
/* Recurse. */
runner_do_unskip(cp, e, drift);
runner_do_unskip(cp, e, 0);
#if 0
dx_max = max(dx_max, cp->dx_max);
h_max = max(h_max, cp->h_max);
......@@ -892,7 +900,7 @@ void runner_do_unskip_mapper(void *map_data, int num_elements,
#ifdef WITH_MPI
if (c != NULL) runner_do_unskip(c, e, (c->nodeID == e->nodeID));
#else
if (c != NULL) runner_do_unskip(c, e, 1);
if (c != NULL) runner_do_unskip(c, e, 0);
#endif
}
}
......
src/runner_doiact.h
View file @ 0b6ee946
......@@ -721,10 +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;
#ifdef SWIFT_DEBUG_CHECKS
cell_is_drifted(ci, e);
cell_is_drifted(cj, e);
#endif
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};
......@@ -764,6 +762,8 @@ 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,6 +774,7 @@ 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];
......@@ -826,6 +827,8 @@ 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;
......@@ -836,6 +839,8 @@ 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];
......@@ -918,11 +923,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;
#ifdef SWIFT_DEBUG_CHECKS
cell_is_drifted(ci, e);
cell_is_drifted(cj, e);
#endif
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);
......@@ -990,6 +993,8 @@ 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;
......@@ -1004,6 +1009,8 @@ 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];
......@@ -1055,6 +1062,8 @@ 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];
......@@ -1134,6 +1143,8 @@ 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;
......@@ -1148,6 +1159,8 @@ 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];
......@@ -1198,6 +1211,8 @@ 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];
......@@ -1313,9 +1328,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
if (!cell_is_active(c, e)) return;
#ifdef SWIFT_DEBUG_CHECKS
cell_is_drifted(c, e);
#endif
if(!cell_is_drifted(c, e)) cell_drift(c, e);
struct part *restrict parts = c->parts;
const int count = c->count;
......@@ -1344,6 +1357,8 @@ 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)) {
......@@ -1354,6 +1369,8 @@ 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];
......@@ -1408,6 +1425,8 @@ 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];
......@@ -1548,9 +1567,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
if (!cell_is_active(c, e)) return;
#ifdef SWIFT_DEBUG_CHECKS
cell_is_drifted(c, e);
#endif
if(!cell_is_drifted(c, e)) error("Cell is not drifted");
struct part *restrict parts = c->parts;
const int count = c->count;
......@@ -1579,6 +1596,8 @@ 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)) {
......@@ -1589,6 +1608,8 @@ 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];
......@@ -1643,6 +1664,8 @@ 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];
......
src/scheduler.c
View file @ 0b6ee946
......@@ -215,7 +215,8 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
/* Get the sort ID, use space_getsid and not t->flags
to make sure we get ci and cj swapped if needed. */
double shift[3];
int sid = space_getsid(s->space, &ci, &cj, shift);
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 &&
......@@ -584,6 +585,7 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
}
/* Otherwise, break it up if it is too large? */
} else if (scheduler_doforcesplit && ci->split && cj->split &&
(ci->count > space_maxsize / cj->count)) {
......@@ -605,6 +607,7 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
}
/* Otherwise, if not spilt, stitch-up the sorting and drift. */
} else {
/* Create the sort for ci. */
......@@ -619,10 +622,21 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
/* Create the drift for ci. */
if (ci->drift == NULL) {
ci->drift = scheduler_addtask(s, task_type_drift, task_subtype_none,
0, 0, ci, NULL, 0);
scheduler_addunlock(s, ci->drift, ci->sorts);
// 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. */
......@@ -637,10 +651,20 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
/* Create the drift for cj. */
if (cj->drift == NULL) {
cj->drift = scheduler_addtask(s, task_type_drift, task_subtype_none,
0, 0, cj, NULL, 0);
scheduler_addunlock(s, cj->drift, cj->sorts);
// 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);
}
......
src/space.c
View file @ 0b6ee946
......@@ -168,6 +168,32 @@ 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.
*
......@@ -346,7 +372,7 @@ void space_regrid(struct space *s, int verbose) {
c->depth = 0;
c->count = 0;
c->gcount = 0;
c->super = c;
// c->super = c;
c->ti_old = ti_current;
lock_init(&c->lock);
}
......@@ -419,6 +445,7 @@ void space_regrid(struct space *s, int verbose) {
s->cells_top[k].cooling = NULL;
s->cells_top[k].sourceterms = NULL;
s->cells_top[k].super = &s->cells_top[k];
s->cells_top[k].ti_old = 0;
#if WITH_MPI
s->cells_top[k].recv_xv = NULL;
s->cells_top[k].recv_rho = NULL;
......@@ -1464,7 +1491,6 @@ void space_split_recursive(struct space *s, struct cell *c, int *buff) {
struct part *parts = c->parts;
struct gpart *gparts = c->gparts;
struct xpart *xparts = c->xparts;
struct engine *e = s->e;
/* If the buff is NULL, allocate it, and remember to free it. */
const int allocate_buffer = (buff == NULL);
......@@ -1494,7 +1520,7 @@ void space_split_recursive(struct space *s, struct cell *c, int *buff) {
temp = space_getcell(s);
temp->count = 0;
temp->gcount = 0;
temp->ti_old = e->ti_current;
temp->ti_old = c->ti_old;
temp->loc[0] = c->loc[0];
temp->loc[1] = c->loc[1];
temp->loc[2] = c->loc[2];
......
src/space.h
View file @ 0b6ee946
......@@ -150,6 +150,7 @@ 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,
......
src/task.c
View file @ 0b6ee946
......@@ -265,7 +265,7 @@ void task_unlock(struct task *t) {
cell_unlocktree(ci);
cell_gunlocktree(ci);
break;
case task_type_sort:
cell_unlocktree(ci);
break;
......@@ -337,8 +337,8 @@ int task_lock(struct task *t) {
if (ci->hold || ci->ghold) return 0;
if (cell_locktree(ci) != 0) return 0;
if (cell_glocktree(ci) != 0) {
cell_unlocktree(ci);
return 0;
cell_unlocktree(ci);
return 0;
}
break;
......
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