diff --git a/examples/main.c b/examples/main.c
index 2c1087db4a1508828cd076de95e5b79811058b5e..ec513723c08322f68395085e122cda9167966444 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -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)
diff --git a/src/engine.c b/src/engine.c
index 53aa744f7664a6dd46962727fcab8290cb36cf16..5819e820dcc6a592134b7e248a2468fa2c5da29c 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -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.
*
diff --git a/src/engine.h b/src/engine.h
index afc2f9ad157b5895c84811d6b61a9e7288131b7b..7c8f03d3813951b4a412545b1f215aaa91bf6dc5 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -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);
diff --git a/src/runner.c b/src/runner.c
index 9544450cd1de0446b125a3e01a880cf67bb8eb04..a7214cad5d9002a6cbe7e6860746ccc817f020f5 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -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) {
diff --git a/src/space.c b/src/space.c
index 2c0e6a7ec36116b28f29c1f728e50b5149eaf651..3ec12b1e2a40b546c8f1c27ea587f5a0ee31c84c 100644
--- a/src/space.c
+++ b/src/space.c
@@ -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.
*
diff --git a/src/space.h b/src/space.h
index cace0a6d6d69e35990309934586c63d0f402f814..a8f51804408566de6f4449b78f6ff2e471dac475 100644
--- a/src/space.h
+++ b/src/space.h
@@ -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);