Skip to content
Snippets Groups Projects
Commit eeb72b2e authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

identified that the bug is in the links between cells in doghost

parent 42c9b1da
No related branches found
No related tags found
2 merge requests!136Master,!79First version of the multiple time-stepping
Hide whitespace changes
Inline Side-by-side
examples/main.c
+ 1
− 2
View file @ eeb72b2e
......@@ -378,7 +378,6 @@ int main(int argc, char *argv[]) {
if (myrank == 0)
printf("# Step Time time-step CPU Wall-clock time [ms]\n");
return 0;
/* Let loose a runner on the space. */
for (j = 0; j < runs && e.time < clock; j++) {
......@@ -397,7 +396,7 @@ int main(int argc, char *argv[]) {
/* Take a step. */
engine_step(&e);
if (j == 0) break;
if (j == 1) break;
if (with_outputs && j % 100 == 0) {
......
src/const.h
+ 1
− 0
View file @ eeb72b2e
......@@ -48,6 +48,7 @@
#define const_eta_kernel \
1.2349f /* Corresponds to 48 ngbs with the cubic spline kernel */
#define const_delta_nwneigh 1.f
#define const_smoothing_max_iter 3
#define CUBIC_SPLINE_KERNEL
/* Gravity stuff. */
......
src/engine.c
+ 1
− 1
View file @ eeb72b2e
......@@ -1252,7 +1252,7 @@ int engine_marktasks(struct engine *e) {
struct scheduler *s = &e->sched;
int k, nr_tasks = s->nr_tasks, *ind = s->tasks_ind;
struct task *t, *tasks = s->tasks;
float t_end = e->time;
//float t_end = e->time;
struct cell *ci, *cj;
// ticks tic = getticks();
......
src/runner.c
+ 35
− 15
View file @ eeb72b2e
......@@ -53,8 +53,8 @@
#define PRINT_PART if(p->id==1000) { \
message("t->t_end p->id=%lld p->h=%f p->N_ngb=%f p->rho=%f p->t_beg=%f p->t_end=%f",\
p->id, p->h, p->density.wcount, p->rho, p->t_begin, p->t_end); \
message("p->id=%lld p->h=%f p->N_ngb=%f p->rho=%f p->t_beg=%f p->t_end=%f",\
p->id, p->h, p->density.wcount, p->rho, p->t_begin, p->t_end); \
}
......@@ -556,7 +556,7 @@ void runner_doghost(struct runner *r, struct cell *c) {
struct part *p, *parts = c->parts;
struct xpart *xp, *xparts = c->xparts;
struct cell *finger;
int i, k, redo, count = c->count;
int i, k, redo, count = c->count, num_reruns;
int *pid;
float h, ih, ih2, ih4, h_corr, rho, wcount, rho_dh, wcount_dh, u, fc;
float normDiv_v, normCurl_v;
......@@ -582,7 +582,7 @@ void runner_doghost(struct runner *r, struct cell *c) {
for (k = 0; k < count; k++) pid[k] = k;
/* While there are particles that need to be updated... */
while (count > 0) {
for (num_reruns = 0; count > 0 && num_reruns < const_smoothing_max_iter; num_reruns++) {
// message("count=%d redo=%d", count, redo);
......@@ -613,7 +613,10 @@ void runner_doghost(struct runner *r, struct cell *c) {
wcount_dh =
p->density.wcount_dh * ih * (4.0f / 3.0 * M_PI * kernel_gamma3);
PRINT_PART
PRINT_PART;
//if(p->id==1000)
// message("wcount_dh=%f", wcount_dh);
/* If no derivative, double the smoothing length. */
if (wcount_dh == 0.0f) h_corr = p->h;
......@@ -634,6 +637,8 @@ void runner_doghost(struct runner *r, struct cell *c) {
/* Ok, correct then */
p->h += h_corr;
//message("Not converged: wcount=%f", p->density.wcount);
/* And flag for another round of fun */
pid[redo] = pid[i];
redo += 1;
......@@ -646,7 +651,7 @@ void runner_doghost(struct runner *r, struct cell *c) {
continue;
}
/* We now have a particle whose smoothing length has converged */
/* We now have a particle whose smoothing length has converged */
/* Pre-compute some stuff for the balsara switch. */
normDiv_v = fabs(p->density.div_v / rho * ih4);
......@@ -703,14 +708,18 @@ void runner_doghost(struct runner *r, struct cell *c) {
count = redo;
if (count > 0) {
// message("count=%d", count);
message("count=%d", count);fflush(stdout);
/* Climb up the cell hierarchy. */
for (finger = c; finger != NULL; finger = finger->parent) {
message("aa"); fflush(stdout);
/* Run through this cell's density interactions. */
for (struct link *l = finger->density; l != NULL; l = l->next) {
message("link: %p next: %p", l, l->next); fflush(stdout);
/* Self-interaction? */
if (l->t->type == task_type_self)
runner_doself_subset_density(r, finger, parts, pid, count);
......@@ -718,6 +727,8 @@ void runner_doghost(struct runner *r, struct cell *c) {
/* Otherwise, pair interaction? */
else if (l->t->type == task_type_pair) {
message("pair");
/* Left or right? */
if (l->t->ci == finger)
runner_dopair_subset_density(r, finger, parts, pid, count,
......@@ -731,6 +742,8 @@ void runner_doghost(struct runner *r, struct cell *c) {
/* Otherwise, sub interaction? */
else if (l->t->type == task_type_sub) {
message("sub");
/* Left or right? */
if (l->t->ci == finger)
runner_dosub_subset_density(r, finger, parts, pid, count,
......@@ -740,10 +753,16 @@ void runner_doghost(struct runner *r, struct cell *c) {
l->t->ci, -1, 1);
}
}
error("done");
}
}
}
if (count)
message("Smoothing length failed to converge on %i particles.", count );
#ifdef TIMER_VERBOSE
message("runner %02i: %i parts at depth %i took %.3f ms.", r->id, c->count,
c->depth, ((double)TIMER_TOC(timer_doghost)) / CPU_TPS * 1000);
......@@ -808,13 +827,13 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
u = p->u *= expf(w);
/* Predict smoothing length */
w = p->force.h_dt * ih * dt;
if (fabsf(w) < 0.01f)
h = p->h *=
1.0f +
w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f * w)));
else
h = p->h *= expf(w);
//w = p->force.h_dt * ih * dt;
//if (fabsf(w) < 0.01f)
// h = p->h *=
// 1.0f +
// w * (1.0f + w * (0.5f + w * (1.0f / 6.0f + 1.0f / 24.0f * w)));
//else
// h = p->h *= expf(w);
/* Predict density */
w = -3.0f * p->force.h_dt * ih * dt;
......@@ -1095,6 +1114,7 @@ void *runner_main(void *data) {
/* Different types of tasks... */
switch (t->type) {
case task_type_self:
//message("self");
if (t->subtype == task_subtype_density)
runner_doself1_density(r, ci);
else if (t->subtype == task_subtype_force)
......
src/scheduler.c
+ 21
− 19
View file @ eeb72b2e
......@@ -496,6 +496,8 @@ void scheduler_splittasks(struct scheduler *s) {
/* Otherwise, if not spilt, stitch-up the sorting. */
else {
//message("called");
/* Create the sort for ci. */
// lock_lock( &ci->lock );
if (ci->sorts == NULL)
......@@ -670,13 +672,15 @@ struct task *scheduler_addtask(struct scheduler *s, int type, int subtype,
/* Get the next free task. */
ind = atomic_inc(&s->tasks_next);
/* Overflow? */
if (ind >= s->size) error("Task list overflow.");
/* Get a pointer to the new task. */
t = &s->tasks[ind];
if (t->type == task_type_sort) message("sort!");
/* Copy the data. */
t->type = type;
t->subtype = subtype;
......@@ -904,8 +908,8 @@ void scheduler_start(struct scheduler *s, unsigned int mask) {
struct task *t, *tasks = s->tasks;
// ticks tic;
message("begin");
fflush(stdout);
//message("begin");
//fflush(stdout);
/* Store the mask */
s->mask = mask;
......@@ -921,17 +925,19 @@ void scheduler_start(struct scheduler *s, unsigned int mask) {
if (!((1 << t->type) & s->mask) || t->skip) continue;
for (j = 0; j < t->nr_unlock_tasks; j++) {
atomic_inc(&t->unlock_tasks[j]->wait);
if(t->unlock_tasks[j] == &tasks[9563] ) {
message("task %d %s %s unlocking task %d %s %s\n",
k, taskID_names[t->type], subtaskID_names[t->subtype],
9563, taskID_names[t->unlock_tasks[j]->type], subtaskID_names[t->unlock_tasks[j]->type]);
}
/* if(t->unlock_tasks[j] == &tasks[9563] ) { */
/* message("task %d %s %s unlocking task %d %s %s\n", */
/* k, taskID_names[t->type], subtaskID_names[t->subtype], */
/* 9563, taskID_names[t->unlock_tasks[j]->type], subtaskID_names[t->unlock_tasks[j]->type]); */
/* } */
}
}
for (k = nr_tasks - 1; k >= 0; k--) {
t = &tasks[tid[k]];
//if (t->type == task_type_sort)
// message("%d %s %s %d %d %d\n", k, taskID_names[t->type], subtaskID_names[t->subtype], t->nr_unlock_tasks, t->wait, t->skip);
if (!((1 << t->type) & s->mask) || t->skip) continue;
fprintf(file, "%d %s %s %d %d\n", k, taskID_names[t->type], subtaskID_names[t->subtype], t->nr_unlock_tasks, t->wait);
}
......@@ -941,7 +947,7 @@ void scheduler_start(struct scheduler *s, unsigned int mask) {
fclose(file);
message("All waits set");
//message("All waits set");
fflush(stdout);
/* Don't enqueue link tasks directly. */
......@@ -961,10 +967,10 @@ void scheduler_start(struct scheduler *s, unsigned int mask) {
}
scheduler_dump_queue(s);
message("Done enqueieing");fflush(stdout);
//message("Done enqueieing");fflush(stdout);
// message( "enqueueing tasks took %.3f ms." , (double)( getticks() - tic ) /
// CPU_TPS * 1000 );
......@@ -984,8 +990,10 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
int err;
#endif
//message("Enqueuing a %s", taskID_names[t->type]);
//fflush(stdout);
//if(t->type == task_type_pair) {
// message("Enqueuing a %s", taskID_names[t->type]);
// fflush(stdout);
// }
/* Ignore skipped tasks and tasks not in the mask. */
if (t->skip || ((1 << t->type) & ~(s->mask) && t->type != task_type_link) ||
......@@ -1087,9 +1095,6 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
struct task *scheduler_done(struct scheduler *s, struct task *t) {
int k, res;
struct task *t2, *next = NULL;
struct cell *super = t->ci->super;
/* Release whatever locks this task held. */
if (!t->implicit) task_unlock(t);
......@@ -1133,9 +1138,6 @@ struct task *scheduler_done(struct scheduler *s, struct task *t) {
struct task *scheduler_unlock(struct scheduler *s, struct task *t) {
int k, res;
struct task *t2, *next = NULL;
/* Loop through the dependencies and add them to a queue if
they are ready. */
for (int k = 0; k < t->nr_unlock_tasks; k++) {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment