diff --git a/examples/main.c b/examples/main.c index 055fb0d99bf7107eaa753d0f706866610cd21f31..5789f03c257cf24a2e48d608551b21de09b68764 100644 --- a/examples/main.c +++ b/examples/main.c @@ -396,7 +396,7 @@ int main(int argc, char *argv[]) { /* Take a step. */ engine_step(&e); - if (j == 1) break; + if (j == 2) break; if (with_outputs && j % 100 == 0) { @@ -505,7 +505,7 @@ int main(int argc, char *argv[]) { #endif /* Say goodbye. */ - if (myrank == 0) message("done."); + if (myrank == 0) message("done. Bye."); /* All is calm, all is bright. */ return 0; diff --git a/src/const.h b/src/const.h index accc682f3eab78ef90bdc70709f0725bbdab1e03..82f5b74604ba315efbbe01c922cb30b973e56965 100644 --- a/src/const.h +++ b/src/const.h @@ -48,7 +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 const_smoothing_max_iter 30 #define CUBIC_SPLINE_KERNEL /* Gravity stuff. */ diff --git a/src/runner.c b/src/runner.c index bc9b8c38f414236b1d34716953e025e4430591cd..3b613dde252a5619636ae7c5af0a6cf5d7a10717 100644 --- a/src/runner.c +++ b/src/runner.c @@ -51,13 +51,12 @@ #endif #include "runner_iact_grav.h" - -#define PRINT_PART if(p->id==1000) { \ - 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); \ -} - - +#define PRINT_PART \ + if (p->id == 1000) { \ + 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); \ + } /* Convert cell location to ID. */ #define cell_getid(cdim, i, j, k) \ @@ -197,8 +196,8 @@ void runner_dosort(struct runner *r, struct cell *c, int flags, int clock) { TIMER_TIC - // message("sort!"); - + // message("sort!"); + /* Clean-up the flags, i.e. filter out what's already been sorted. */ flags &= ~c->sorted; if (flags == 0) return; @@ -540,7 +539,6 @@ void runner_doinit(struct runner *r, struct cell *c) { } PRINT_PART; - } } @@ -567,7 +565,7 @@ void runner_doghost(struct runner *r, struct cell *c) { TIMER_TIC; - //message("start"); + // message("start"); /* Recurse? */ if (c->split) { @@ -582,7 +580,8 @@ 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... */ - for (num_reruns = 0; count > 0 && num_reruns < const_smoothing_max_iter; num_reruns++) { + for (num_reruns = 0; count > 0 && num_reruns < const_smoothing_max_iter; + num_reruns++) { // message("count=%d redo=%d", count, redo); @@ -613,11 +612,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; - //if(p->id==1000) - // message("wcount_dh=%f", wcount_dh); - - + 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; @@ -637,8 +635,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); - + // message("Not converged: wcount=%f", p->density.wcount); + /* And flag for another round of fun */ pid[redo] = pid[i]; redo += 1; @@ -651,7 +649,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); @@ -708,16 +706,16 @@ void runner_doghost(struct runner *r, struct cell *c) { count = redo; if (count > 0) { - message("count=%d", count);fflush(stdout); + // message("count=%d", count);fflush(stdout); /* Climb up the cell hierarchy. */ for (finger = c; finger != NULL; finger = finger->parent) { - + /* Run through this cell's density interactions. */ for (struct link *l = finger->link_density; l != NULL; l = l->next) { - - //message("link: %p next: %p", l, l->next); fflush(stdout); - + + // 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); @@ -725,8 +723,8 @@ void runner_doghost(struct runner *r, struct cell *c) { /* Otherwise, pair interaction? */ else if (l->t->type == task_type_pair) { - //message("pair"); - + // message("pair"); + /* Left or right? */ if (l->t->ci == finger) runner_dopair_subset_density(r, finger, parts, pid, count, @@ -740,8 +738,8 @@ void runner_doghost(struct runner *r, struct cell *c) { /* Otherwise, sub interaction? */ else if (l->t->type == task_type_sub) { - //message("sub"); - + // message("sub"); + /* Left or right? */ if (l->t->ci == finger) runner_dosub_subset_density(r, finger, parts, pid, count, @@ -751,16 +749,14 @@ void runner_doghost(struct runner *r, struct cell *c) { l->t->ci, -1, 1); } } - error("done"); + // error("done"); } } } - - if (count) - message("Smoothing length failed to converge on %i particles.", count ); + 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); @@ -801,7 +797,7 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) { ih = 1.0f / h; PRINT_PART; - + /* Drift... */ p->x[0] += xp->v_full[0] * dt; p->x[1] += xp->v_full[1] * dt; @@ -825,12 +821,12 @@ 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) + // 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 + // else // h = p->h *= expf(w); /* Predict density */ @@ -940,8 +936,8 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) { p->t_begin = p->t_end; p->t_end = p->t_begin + dt; - PRINT_PART - + PRINT_PART + /* Kick particles in momentum space */ xp->v_full[0] += p->a[0] * dt; xp->v_full[1] += p->a[1] * dt; @@ -1068,8 +1064,9 @@ void *runner_main(void *data) { struct part *parts; int k, nr_parts; - message("here");fflush(stdout); - + message("here"); + fflush(stdout); + /* Main loop. */ while (1) { @@ -1090,8 +1087,8 @@ void *runner_main(void *data) { t = scheduler_gettask(sched, r->qid, super); TIMER_TOC(timer_gettask); - //message("Got task %p", t->type); fflush(stdout); - + // message("Got task %p", t->type); fflush(stdout); + /* Did I get anything? */ if (t == NULL) break; } @@ -1112,7 +1109,7 @@ void *runner_main(void *data) { /* Different types of tasks... */ switch (t->type) { case task_type_self: - //message("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 14c6f2820f0a364088aba89169c26f9e6824bf15..a5a7c3dd19ab3966408dfb4dae3cfde3cf74e442 100644 --- a/src/scheduler.c +++ b/src/scheduler.c @@ -496,8 +496,8 @@ void scheduler_splittasks(struct scheduler *s) { /* Otherwise, if not spilt, stitch-up the sorting. */ else { - //message("called"); - + // message("called"); + /* Create the sort for ci. */ // lock_lock( &ci->lock ); if (ci->sorts == NULL) @@ -672,7 +672,7 @@ 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."); @@ -680,7 +680,7 @@ struct task *scheduler_addtask(struct scheduler *s, int type, int subtype, t = &s->tasks[ind]; if (t->type == task_type_sort) message("sort!"); - + /* Copy the data. */ t->type = type; t->subtype = subtype; @@ -908,14 +908,14 @@ 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; - FILE* file = fopen("tasks.dat", "w"); - + FILE *file = fopen("tasks.dat", "w"); + /* Run through the tasks and set their waits. */ // tic = getticks(); for (k = nr_tasks - 1; k >= 0; k--) { @@ -927,29 +927,32 @@ void scheduler_start(struct scheduler *s, unsigned int mask) { 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]); */ + /* 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); - } + 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); + } // message( "waiting tasks took %.3f ms." , (double)( getticks() - tic ) / // CPU_TPS * 1000 ); fclose(file); - - //message("All waits set"); + + // message("All waits set"); fflush(stdout); - + /* Don't enqueue link tasks directly. */ s->mask &= ~(1 << task_type_link); @@ -967,11 +970,9 @@ 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 ); } @@ -990,11 +991,11 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) { int err; #endif - //if(t->type == task_type_pair) { + // 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) || atomic_cas(&t->rid, -1, 0) != -1) @@ -1095,7 +1096,6 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) { struct task *scheduler_done(struct scheduler *s, struct task *t) { - /* Release whatever locks this task held. */ if (!t->implicit) task_unlock(t); @@ -1287,48 +1287,43 @@ void scheduler_init(struct scheduler *s, struct space *space, int nr_queues, s->tasks_next = 0; } - /** * @brief Print all the tasks in the queue of the scheduler * * @param s The #scheduler. - */ - void scheduler_dump_queue(struct scheduler *s) { - - int i,j; - FILE* file; - char buffer[256]; - struct queue* q; - struct task* t; - - for( i=0; i<s->nr_queues; ++i) { - - - /* Open file */ - sprintf(buffer, "queue_%d.dat", i); - file = fopen(buffer, "w"); - - /* Get the queue */ - q = &s->queues[i]; - - /* Some general info */ - fprintf(file, "# Queue %d, size=%d, count=%d\n", i, q->size, q->count); - fprintf(file, "# Index type subtype\n"); - - for (j=0; j<q->count; ++j) { - - /* Get the task */ - t = &q->tasks[j]; - - /* And print... */ - fprintf(file, "%d %s %s\n", j , taskID_names[t->type], subtaskID_names[t->subtype]); - - } - - - - /* Be nice and clean */ - fclose(file); - } - - } + */ +void scheduler_dump_queue(struct scheduler *s) { + + int i, j; + FILE *file; + char buffer[256]; + struct queue *q; + struct task *t; + + for (i = 0; i < s->nr_queues; ++i) { + + /* Open file */ + sprintf(buffer, "queue_%d.dat", i); + file = fopen(buffer, "w"); + + /* Get the queue */ + q = &s->queues[i]; + + /* Some general info */ + fprintf(file, "# Queue %d, size=%d, count=%d\n", i, q->size, q->count); + fprintf(file, "# Index type subtype\n"); + + for (j = 0; j < q->count; ++j) { + + /* Get the task */ + t = &q->tasks[j]; + + /* And print... */ + fprintf(file, "%d %s %s\n", j, taskID_names[t->type], + subtaskID_names[t->subtype]); + } + + /* Be nice and clean */ + fclose(file); + } +} diff --git a/src/task.c b/src/task.c index 665434bd413b94eae850b650b0914e110f5ab91c..8a9ef80886f8f0624215bef6347741c14d84a595 100644 --- a/src/task.c +++ b/src/task.c @@ -48,9 +48,8 @@ const char *taskID_names[task_type_count] = { "ghost", "drift", "kick", "send", "recv", "link", "grav_pp", "grav_mm", "grav_up", "grav_down"}; -const char *subtaskID_names[task_type_count] = { - "none", "density", "force", "grav"}; - +const char *subtaskID_names[task_type_count] = {"none", "density", + "force", "grav"}; /** * @brief Unlock the cell held by this task. diff --git a/src/tools.c b/src/tools.c index 13a1df4187698337d7959918b12ac1407d924177..b6a945ffba39cdba3cbbf460e61133adc6794f97 100644 --- a/src/tools.c +++ b/src/tools.c @@ -37,7 +37,7 @@ void factor(int value, int *f1, int *f2) { int j; int i; - j = (int) sqrt(value); + j = (int)sqrt(value); for (i = j; i > 0; i--) { if ((value % i) == 0) { *f1 = i; @@ -110,7 +110,7 @@ void pairs_n2(double *dim, struct part *__restrict__ parts, int N, /* Dump the result. */ printf("pairs_n2: avg. density per part is %.3f (nr. pairs %.3f).\n", - rho / N + 32.0 / 3, ((double) count) / N); + rho / N + 32.0 / 3, ((double)count) / N); printf("pairs_n2: densities are in [ %e , %e ].\n", rho_min / N + 32.0 / 3, rho_max / N + 32.0 / 3); /* printf( "pairs_n2: maximum ratio between parts %i [%e,%e,%e] and %i @@ -183,7 +183,7 @@ void pairs_single_grav(double *dim, long long int pid, // int mj, mk; // double maxratio = 1.0; double r2, dx[3]; - float fdx[3], a[3] = { 0.0, 0.0, 0.0 }, aabs[3] = { 0.0, 0.0, 0.0 }; + float fdx[3], a[3] = {0.0, 0.0, 0.0}, aabs[3] = {0.0, 0.0, 0.0}; struct gpart pi, pj; // double ih = 12.0/6.25; @@ -239,7 +239,7 @@ void pairs_single_grav(double *dim, long long int pid, void density_dump(int N) { int k; - float r2[4] = { 0.0f, 0.0f, 0.0f, 0.0f }, hi[4], hj[4]; + float r2[4] = {0.0f, 0.0f, 0.0f, 0.0f}, hi[4], hj[4]; struct part *pi[4], *pj[4], Pi[4], Pj[4]; /* Init the interaction parameters. */ @@ -260,7 +260,7 @@ void density_dump(int N) { r2[3] = r2[2]; r2[2] = r2[1]; r2[1] = r2[0]; - r2[0] = ((float) k) / N; + r2[0] = ((float)k) / N; Pi[0].density.wcount = 0; Pj[0].density.wcount = 0; runner_iact_density(r2[0], NULL, hi[0], hj[0], &Pi[0], &Pj[0]); @@ -305,8 +305,7 @@ void engine_single_density(double *dim, long long int pid, p.density.wcount = 0.0f; p.density.wcount_dh = 0.0f; p.density.div_v = 0.0; - for (k = 0; k < 3; k++) - p.density.curl_v[k] = 0.0; + for (k = 0; k < 3; k++) p.density.curl_v[k] = 0.0; /* Loop over all particle pairs (force). */ for (k = 0; k < N; k++) { @@ -335,7 +334,6 @@ void engine_single_density(double *dim, long long int pid, message("part %lli (h=%e) has wcount=%e, rho=%e, rho_dh=%e.", p.id, p.h, p.density.wcount, p.rho, p.rho_dh); fflush(stdout); - } void engine_single_force(double *dim, long long int pid, @@ -388,8 +386,7 @@ void engine_single_force(double *dim, long long int pid, } /* Dump the result. */ - message( "part %lli (h=%e) has a=[%.3e,%.3e,%.3e], udt=%e." , p.id , - p.h , p.a[0] , p.a[1] , p.a[2] , p.force.u_dt ); + message("part %lli (h=%e) has a=[%.3e,%.3e,%.3e], udt=%e.", p.id, p.h, p.a[0], + p.a[1], p.a[2], p.force.u_dt); fflush(stdout); - }