diff --git a/examples/main.c b/examples/main.c
index 717ae74b27d1b997418dcf61a860a844bd049a3f..055fb0d99bf7107eaa753d0f706866610cd21f31 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -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) {
diff --git a/src/const.h b/src/const.h
index d62673332c2d19872cad1d9c62ca6eff7f0e013a..accc682f3eab78ef90bdc70709f0725bbdab1e03 100644
--- a/src/const.h
+++ b/src/const.h
@@ -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. */
diff --git a/src/engine.c b/src/engine.c
index 666c7d0aa6a010f0cb54d629e39b29c4e7b53311..1018552a1fa2f7a122be4994d634811a578b9cb3 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -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();
diff --git a/src/runner.c b/src/runner.c
index d0b2b2f85720448c56901524131788d163eeda9a..83a9508dd427d8f3dba161cbaca4fc2246bcdf15 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -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)
diff --git a/src/scheduler.c b/src/scheduler.c
index 2b7bcd8cf0a74f05b13abb8b228f276cbd1d6795..14c6f2820f0a364088aba89169c26f9e6824bf15 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -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++) {