Skip to content
GitLab
Commit a957d0a9 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

We now have a proper initialisation step.

parent 78a1b843
Hide whitespace changes
Inline Side-by-side
examples/main.c
View file @ a957d0a9
......@@ -85,8 +85,8 @@ int main(int argc, char *argv[]) {
FILE *file_thread;
int with_outputs = 1;
/* Choke on FP-exceptions. */
// feenableexcept( FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW );
/* Choke on FP-exceptions. */
//feenableexcept( FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW );
#ifdef WITH_MPI
/* Start by initializing MPI. */
......@@ -207,6 +207,7 @@ int main(int argc, char *argv[]) {
/* How large are the parts? */
if (myrank == 0) {
message("sizeof(struct part) is %li bytes.", (long int)sizeof(struct part));
message("sizeof(struct xpart) is %li bytes.", (long int)sizeof(struct xpart));
message("sizeof(struct gpart) is %li bytes.",
(long int)sizeof(struct gpart));
}
......@@ -370,6 +371,9 @@ int main(int argc, char *argv[]) {
fflush(stdout);
}
/* Initialise the particles */
engine_init_particles(&e);
/* Legend */
if (myrank == 0)
printf("# Step Time time-step CPU Wall-clock time [ms]\n");
......@@ -390,6 +394,8 @@ int main(int argc, char *argv[]) {
/* Take a step. */
engine_step(&e);
break;
if (with_outputs && j % 100 == 0) {
#if defined(WITH_MPI)
......
src/engine.c
View file @ a957d0a9
......@@ -93,12 +93,16 @@ void engine_mkghosts(struct engine *e, struct cell *c, struct cell *super) {
int k;
struct scheduler *s = &e->sched;
// message("in here");
/* Am I the super-cell? */
if (super == NULL && c->nr_tasks > 0) {
/* Remember me. */
super = c;
//message("Adding tasks");
/* Local tasks only... */
if (c->nodeID == e->nodeID) {
......@@ -1148,6 +1152,7 @@ void engine_maketasks(struct engine *e) {
/* Self-interaction? */
if (t->type == task_type_self && t->subtype == task_subtype_density) {
scheduler_addunlock(sched, t->ci->super->init, t);
scheduler_addunlock(sched, t, t->ci->super->ghost);
t2 = scheduler_addtask(sched, task_type_self, task_subtype_force, 0, 0,
t->ci, NULL, 0);
......@@ -1162,11 +1167,13 @@ void engine_maketasks(struct engine *e) {
t2 = scheduler_addtask(sched, task_type_pair, task_subtype_force, 0, 0,
t->ci, t->cj, 0);
if (t->ci->nodeID == nodeID) {
scheduler_addunlock(sched, t->ci->super->init, t);
scheduler_addunlock(sched, t, t->ci->super->ghost);
scheduler_addunlock(sched, t->ci->super->ghost, t2);
scheduler_addunlock(sched, t2, t->ci->super->kick);
}
if (t->cj->nodeID == nodeID && t->ci->super != t->cj->super) {
scheduler_addunlock(sched, t->cj->super->init, t);
scheduler_addunlock(sched, t, t->cj->super->ghost);
scheduler_addunlock(sched, t->cj->super->ghost, t2);
scheduler_addunlock(sched, t2, t->cj->super->kick);
......@@ -1182,12 +1189,14 @@ void engine_maketasks(struct engine *e) {
t2 = scheduler_addtask(sched, task_type_sub, task_subtype_force, t->flags,
0, t->ci, t->cj, 0);
if (t->ci->nodeID == nodeID) {
scheduler_addunlock(sched, t->ci->super->init, t);
scheduler_addunlock(sched, t, t->ci->super->ghost);
scheduler_addunlock(sched, t->ci->super->ghost, t2);
scheduler_addunlock(sched, t2, t->ci->super->kick);
}
if (t->cj != NULL && t->cj->nodeID == nodeID &&
t->ci->super != t->cj->super) {
scheduler_addunlock(sched, t->cj->super->init, t);
scheduler_addunlock(sched, t, t->cj->super->ghost);
scheduler_addunlock(sched, t->cj->super->ghost, t2);
scheduler_addunlock(sched, t2, t->cj->super->kick);
......@@ -1252,116 +1261,119 @@ int engine_marktasks(struct engine *e) {
struct cell *ci, *cj;
// ticks tic = getticks();
/* Much less to do here if we're on a fixed time-step. */
if (!(e->policy & engine_policy_multistep)) {
/* /\* Much less to do here if we're on a fixed time-step. *\/ */
/* if (!(e->policy & engine_policy_multistep)) { */
/* Run through the tasks and mark as skip or not. */
for (k = 0; k < nr_tasks; k++) {
/* /\* Run through the tasks and mark as skip or not. *\/ */
/* for (k = 0; k < nr_tasks; k++) { */
/* Get a handle on the kth task. */
t = &tasks[ind[k]];
/* /\* Get a handle on the kth task. *\/ */
/* t = &tasks[ind[k]]; */
/* Pair? */
if (t->type == task_type_pair ||
(t->type == task_type_sub && t->cj != NULL)) {
/* /\* Pair? *\/ */
/* if (t->type == task_type_pair || */
/* (t->type == task_type_sub && t->cj != NULL)) { */
/* Local pointers. */
ci = t->ci;
cj = t->cj;
/* /\* Local pointers. *\/ */
/* ci = t->ci; */
/* cj = t->cj; */
/* Too much particle movement? */
if (t->tight &&
(fmaxf(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin ||
ci->dx_max > space_maxreldx * ci->h_max ||
cj->dx_max > space_maxreldx * cj->h_max))
return 1;
/* /\* Too much particle movement? *\/ */
/* if (t->tight && */
/* (fmaxf(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin || */
/* ci->dx_max > space_maxreldx * ci->h_max || */
/* cj->dx_max > space_maxreldx * cj->h_max)) */
/* return 1; */
}
/* } */
/* Sort? */
else if (t->type == task_type_sort) {
/* /\* Sort? *\/ */
/* else if (t->type == task_type_sort) { */
/* If all the sorts have been done, make this task implicit. */
if (!(t->flags & (t->flags ^ t->ci->sorted))) t->implicit = 1;
}
}
/* /\* If all the sorts have been done, make this task implicit. *\/ */
/* if (!(t->flags & (t->flags ^ t->ci->sorted))) t->implicit = 1; */
/* } */
/* } */
} else {
/* } else { */
/* Run through the tasks and mark as skip or not. */
for (k = 0; k < nr_tasks; k++) {
/* Run through the tasks and mark as skip or not. */
for (k = 0; k < nr_tasks; k++) {
/* Get a handle on the kth task. */
t = &tasks[ind[k]];
/* Get a handle on the kth task. */
t = &tasks[ind[k]];
/* Sort-task? Note that due to the task ranking, the sorts
will all come before the pairs. */
if (t->type == task_type_sort) {
/* Sort-task? Note that due to the task ranking, the sorts
will all come before the pairs. */
if (t->type == task_type_sort) {
/* Re-set the flags. */
t->flags = 0;
t->skip = 1;
/* Re-set the flags. */
t->flags = 0;
t->skip = 1;
}
/* Single-cell task? */
else if (t->type == task_type_self || t->type == task_type_ghost ||
(t->type == task_type_sub && t->cj == NULL)) {
}
/* Set this task's skip. */
t->skip = (t->ci->t_end_min > t_end);
/* Single-cell task? */
else if (t->type == task_type_self || t->type == task_type_ghost ||
(t->type == task_type_sub && t->cj == NULL)) {
}
/* Set this task's skip. */
//t->skip = (t->ci->t_end_min >= t_end);
/* Pair? */
else if (t->type == task_type_pair ||
(t->type == task_type_sub && t->cj != NULL)) {
/* Local pointers. */
ci = t->ci;
cj = t->cj;
/* Set this task's skip. */
t->skip = (ci->t_end_min > t_end && cj->t_end_min > t_end);
/* Too much particle movement? */
if (t->tight &&
(fmaxf(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin ||
ci->dx_max > space_maxreldx * ci->h_max ||
cj->dx_max > space_maxreldx * cj->h_max))
return 1;
/* Set the sort flags. */
if (!t->skip && t->type == task_type_pair) {
if (!(ci->sorted & (1 << t->flags))) {
ci->sorts->flags |= (1 << t->flags);
ci->sorts->skip = 0;
}
if (!(cj->sorted & (1 << t->flags))) {
cj->sorts->flags |= (1 << t->flags);
cj->sorts->skip = 0;
}
}
}
/* Pair? */
else if (t->type == task_type_pair ||
(t->type == task_type_sub && t->cj != NULL)) {
/* Local pointers. */
ci = t->ci;
cj = t->cj;
/* Set this task's skip. */
//t->skip = (ci->t_end_min >= t_end && cj->t_end_min >= t_end);
/* Too much particle movement? */
if (t->tight &&
(fmaxf(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin ||
ci->dx_max > space_maxreldx * ci->h_max ||
cj->dx_max > space_maxreldx * cj->h_max))
return 1;
/* Set the sort flags. */
if (!t->skip && t->type == task_type_pair) {
if (!(ci->sorted & (1 << t->flags))) {
ci->sorts->flags |= (1 << t->flags);
ci->sorts->skip = 0;
}
if (!(cj->sorted & (1 << t->flags))) {
cj->sorts->flags |= (1 << t->flags);
cj->sorts->skip = 0;
}
}
/* Kick? */
else if (t->type == task_type_kick)
t->skip = 0;
}
/* Drift? */
else if (t->type == task_type_drift)
t->skip = 0;
/* Kick? */
else if (t->type == task_type_kick)
t->skip = 0;
/* Init? */
else if (t->type == task_type_init)
t->skip = 0;
/* Drift? */
else if (t->type == task_type_drift)
t->skip = 0;
/* None? */
else if (t->type == task_type_none)
t->skip = 1;
}
/* Init? */
else if (t->type == task_type_init)
t->skip = 0;
/* None? */
else if (t->type == task_type_none)
t->skip = 1;
}
//}
// message( "took %.3f ms." , (double)(getticks() - tic)/CPU_TPS*1000 );
......@@ -1742,6 +1754,78 @@ void hassorted(struct cell *c) {
if (c->progeny[k] != NULL) hassorted(c->progeny[k]);
}
/**
* @brief Initialises the particles and set them in a state ready to move forward in time.
*
* @param e The #engine
*/
void engine_init_particles(struct engine *e) {
struct space *s = e->s;
//engine_repartition(e);
engine_prepare(e);
engine_print(e);
//engine_maketasks(e);
engine_print(e);
engine_marktasks(e);
engine_print(e);
fflush(stdout);
message("Engine prepared");
/* Make sure all particles are ready to go */
void initParts(struct part *p, struct xpart *xp, struct cell *c) {
p->t_end = 0.;
p->t_begin = 0.;
p->rho = -1.;
p->h = 1.12349f / 20;
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
c->t_end_min = 0.;
}
message("Initialising particles");
space_map_parts_xparts(s, initParts);
/* Now everybody should have sensible smoothing length */
void printParts(struct part *p, struct xpart *xp, struct cell *c) {
if( p->id == 1000)
message("id=%lld h=%f rho=%f", p->id, p->h, p->rho);
}
//space_map_parts_xparts(s, printParts);
void printCells(struct part *p, struct xpart *xp, struct cell *c) {
if(c->super != NULL && 0)
message("c->t_end_min=%f c->t_end_max=%f c->super=%p sort=%p ghost=%p kick=%p", c->t_end_min, c->t_end_max, c->super, c->sorts, c->ghost, c->kick);
}
space_map_parts_xparts(s, printCells);
/* Now do a density calculation */
TIMER_TIC;
engine_launch(e, e->nr_threads,
(1 << task_type_sort) | (1 << task_type_self) |
(1 << task_type_pair) | (1 << task_type_sub) |
(1 << task_type_init) | (1 << task_type_ghost) |
(1 << task_type_send) | (1 << task_type_recv) |
(1 << task_type_link));
TIMER_TOC(timer_runners);
space_map_parts_xparts(s, printParts);
}
/**
* @brief Let the #engine loose to compute the forces.
*
......
src/engine.h
View file @ a957d0a9
......@@ -138,6 +138,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
float timeBegin, float timeEnd);
void engine_prepare(struct engine *e);
void engine_print(struct engine *e);
void engine_init_particles(struct engine *e);
void engine_step(struct engine *e);
void engine_maketasks(struct engine *e);
void engine_split(struct engine *e, int *grid);
......
src/runner.c
View file @ a957d0a9
......@@ -512,6 +512,8 @@ void runner_doinit(struct runner *r, struct cell *c) {
return;
}
// message("in here count=%d", count);
/* Loop over the parts in this cell. */
for (i = 0; i < count; i++) {
......@@ -528,6 +530,12 @@ void runner_doinit(struct runner *r, struct cell *c) {
p->density.div_v = 0.0;
for (k = 0; k < 3; k++) p->density.curl_v[k] = 0.0;
}
if(p->id==1000) {
message("p->id=%lld p->h=%f p->rho=%f", p->id, p->h, p->rho);
}
}
}
......@@ -554,6 +562,8 @@ void runner_doghost(struct runner *r, struct cell *c) {
TIMER_TIC
//message("start");
/* Recurse? */
if (c->split) {
for (k = 0; k < 8; k++)
......@@ -561,6 +571,8 @@ void runner_doghost(struct runner *r, struct cell *c) {
return;
}
//message("done recursing");
/* Init the IDs that have to be updated. */
if ((pid = (int *)alloca(sizeof(int) * count)) == NULL)
error("Call to alloca failed.");
......@@ -569,9 +581,11 @@ void runner_doghost(struct runner *r, struct cell *c) {
/* While there are particles that need to be updated... */
while (count > 0) {
//message("count=%d redo=%d", count, redo);
/* Reset the redo-count. */
redo = 0;
/* Loop over the parts in this cell. */
for (i = 0; i < count; i++) {
......@@ -596,6 +610,12 @@ void runner_doghost(struct runner *r, struct cell *c) {
wcount_dh =
p->density.wcount_dh * ih * (4.0f / 3.0 * M_PI * kernel_gamma3);
if(p->id==1000) {
message("p->id=%lld p->h=%f p->N_ngb=%f p->rho=%f", p->id, p->h, wcount, p->rho);
//error("");
}
/* If no derivative, double the smoothing length. */
if (wcount_dh == 0.0f) h_corr = p->h;
......@@ -612,8 +632,8 @@ void runner_doghost(struct runner *r, struct cell *c) {
if (wcount > kernel_nwneigh + const_delta_nwneigh ||
wcount < kernel_nwneigh - const_delta_nwneigh) {
h = p->h += h_corr;
p->h += (p->density.wcount + kernel_root - kernel_nwneigh) /
p->density.wcount_dh;
//p->h += (p->density.wcount + kernel_root - kernel_nwneigh) /
// p->density.wcount_dh;
pid[redo] = pid[i];
redo += 1;
p->density.wcount = 0.0;
......@@ -672,9 +692,11 @@ void runner_doghost(struct runner *r, struct cell *c) {
p->force.u_dt = 0.0f;
p->force.h_dt = 0.0f;
p->force.v_sig = 0.0f;
}
}
/* We now need to treat the particles whose smoothing length had not
* converged again */
......@@ -682,8 +704,8 @@ void runner_doghost(struct runner *r, struct cell *c) {
count = redo;
if (count > 0) {
// error( "Bad smoothing length, fixing this isn't implemented yet." );
//message("count=%d", count);
/* Climb up the cell hierarchy. */
for (finger = c; finger != NULL; finger = finger->parent) {
......
src/space.c
View file @ a957d0a9
......@@ -824,6 +824,48 @@ void space_map_parts(struct space *s,
rec_map_parts(&s->cells[cid], fun, data);
}
/**
* @brief Map a function to all particles in a cell recursively.
*
* @param c The #cell we are working in.
* @param fun Function pointer to apply on the cells.
* @param data Data passed to the function fun.
*/
static void rec_map_parts_xparts(struct cell *c,
void (*fun)(struct part *p, struct xpart *xp, struct cell *c)) {
int k;
/* No progeny? */
if (!c->split)
for (k = 0; k < c->count; k++) fun(&c->parts[k], &c->xparts[k], c);
/* Otherwise, recurse. */
else
for (k = 0; k < 8; k++)
if (c->progeny[k] != NULL) rec_map_parts_xparts(c->progeny[k], fun);
}
/**
* @brief Map a function to all particles (#part and #xpart) in a space.
*
* @param s The #space we are working in.
* @param fun Function pointer to apply on the particles in the cells.
*/
void space_map_parts_xparts(struct space *s,
void (*fun)(struct part *p, struct xpart *xp, struct cell *c)) {
int cid = 0;
/* Call the recursive function on all higher-level cells. */
for (cid = 0; cid < s->nr_cells; cid++)
rec_map_parts_xparts(&s->cells[cid], fun);
}
/**
* @brief Map a function to all particles in a cell recursively.
*
......
src/space.h
View file @ a957d0a9
......@@ -125,6 +125,8 @@ void space_map_cells_pre(struct space *s, int full,
void space_map_parts(struct space *s,
void (*fun)(struct part *p, struct cell *c, void *data),
void *data);
void space_map_parts_xparts(struct space *s,
void (*fun)(struct part *p, struct xpart *xp, struct cell *c));
void space_map_cells_post(struct space *s, int full,
void (*fun)(struct cell *c, void *data), void *data);
void space_rebuild(struct space *s, double h_max, int verbose);
......
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