diff --git a/src/engine.c b/src/engine.c
index 5819e820dcc6a592134b7e248a2468fa2c5da29c..589b5f9ecb9da50de2cf3945fee96af6afcbbf2f 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -49,11 +49,6 @@
#include "error.h"
#include "timers.h"
-#ifdef LEGACY_GADGET2_SPH
-#include "runner_iact_legacy.h"
-#else
-#include "runner_iact.h"
-#endif
/* Convert cell location to ID. */
#define cell_getid(cdim, i, j, k) \
@@ -1728,6 +1723,11 @@ void engine_launch(struct engine *e, int nr_runners, unsigned int mask) {
if (pthread_cond_broadcast(&e->barrier_cond) != 0)
error("Failed to broadcast barrier open condition.");
+ /* Print out what we do */
+ printf("\n");
+ task_print_mask(mask);
+ printf("\n");
+
/* Load the tasks. */
pthread_mutex_unlock(&e->barrier_mutex);
scheduler_start(&e->sched, mask);
@@ -1745,14 +1745,14 @@ void engine_launch(struct engine *e, int nr_runners, unsigned int mask) {
error("Error while waiting for barrier.");
}
-void hassorted(struct cell *c) {
+/* void hassorted(struct cell *c) { */
- if (c->sorted) error("Suprious sorted flags.");
+/* if (c->sorted) error("Suprious sorted flags."); */
- if (c->split)
- for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL) hassorted(c->progeny[k]);
-}
+/* if (c->split) */
+/* for (int k = 0; k < 8; k++) */
+/* 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.
@@ -1782,10 +1782,9 @@ void engine_init_particles(struct engine *e) {
/* 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->t_end = 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];
@@ -1799,7 +1798,7 @@ void engine_init_particles(struct engine *e) {
/* 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);
+ message("id=%lld h=%f rho=%f t_begin=%f t_end=%f", p->id, p->h, p->rho, p->t_begin, p->t_end);
}
//space_map_parts_xparts(s, printParts);
@@ -1807,7 +1806,7 @@ void engine_init_particles(struct engine *e) {
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);
+ //space_map_parts_xparts(s, printCells);
/* Now do a density calculation */
@@ -1821,9 +1820,12 @@ void engine_init_particles(struct engine *e) {
TIMER_TOC(timer_runners);
- space_map_parts_xparts(s, printParts);
+ //space_map_parts_xparts(s, printParts);
-
+ printf("\n\n");
+
+ /* Ready to go */
+ e->step = -1;
}
/**
@@ -1844,28 +1846,7 @@ void engine_step(struct engine *e) {
TIMER_TIC2;
- message("Begin step: %d", e->step);
-
- /* Re-distribute the particles amongst the nodes? */
- if (e->forcerepart) engine_repartition(e);
-
- /* Prepare the space. */
- engine_prepare(e);
-
- engine_print(e);
-
- /* Send off the runners. */
- 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_kick) | (1 << task_type_send) |
- (1 << task_type_recv) | (1 << task_type_link));
-
- TIMER_TOC(timer_runners);
-
- /* Collect the cell data from the kick. */
+ /* Collect the cell data. */
for (k = 0; k < s->nr_cells; k++)
if (s->cells[k].nodeID == e->nodeID) {
c = &s->cells[k];
@@ -1905,24 +1886,49 @@ void engine_step(struct engine *e) {
count = in[0];
ekin = in[1];
epot = in[2];
-/* int nr_parts;
-if ( MPI_Allreduce( &s->nr_parts , &nr_parts , 1 , MPI_INT , MPI_SUM ,
-MPI_COMM_WORLD ) != MPI_SUCCESS )
- error( "Failed to aggregate particle count." );
-if ( e->nodeID == 0 )
- message( "nr_parts=%i." , nr_parts ); */
#endif
+
+ message("t_end_min=%f t_end_max=%f", t_end_min, t_end_max);
+
+
/* Move forward in time */
e->timeOld = e->time;
e->time = t_end_min;
e->step += 1;
+ message("Step: %d e->time=%f", e->step, e->time);
+
+
+ /* Drift everybody */
+ engine_launch(e, e->nr_threads, 1<<task_type_drift);
+
+
+ /* Re-distribute the particles amongst the nodes? */
+ if (e->forcerepart) engine_repartition(e);
+
+ /* Prepare the space. */
+ engine_prepare(e);
+
+ engine_maketasks(e);
+
+ engine_marktasks(e);
+
+ engine_print(e);
+ message("Go !");
+ fflush(stdout);
+
+ /* Send off the runners. */
+ 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_kick) | (1 << task_type_send) |
+ (1 << task_type_recv) | (1 << task_type_link));
+
+ TIMER_TOC(timer_runners);
- message("e->time=%f", e->time);
- message("Ready to drift");
- /* Drift all particles */
- engine_launch(e, e->nr_threads, (1 << task_type_drift));
TIMER_TOC2(timer_step);
}
diff --git a/src/runner.c b/src/runner.c
index a7214cad5d9002a6cbe7e6860746ccc817f020f5..c4837b63b64669e7a558cf75cdb80ae6684f935a 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -51,6 +51,14 @@
#endif
#include "runner_iact_grav.h"
+
+#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); \
+}
+
+
+
/* Convert cell location to ID. */
#define cell_getid(cdim, i, j, k) \
((int)(k) + (cdim)[2] * ((int)(j) + (cdim)[1] * (int)(i)))
@@ -189,6 +197,8 @@ void runner_dosort(struct runner *r, struct cell *c, int flags, int clock) {
TIMER_TIC
+ message("sort!");
+
/* Clean-up the flags, i.e. filter out what's already been sorted. */
flags &= ~c->sorted;
if (flags == 0) return;
@@ -512,8 +522,6 @@ 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++) {
@@ -531,10 +539,7 @@ void runner_doinit(struct runner *r, struct cell *c) {
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);
- }
+ PRINT_PART;
}
}
@@ -560,10 +565,10 @@ void runner_doghost(struct runner *r, struct cell *c) {
#endif
float t_end = r->e->time;
- TIMER_TIC
+ TIMER_TIC;
+
+ //message("start");
- //message("start");
-
/* Recurse? */
if (c->split) {
for (k = 0; k < 8; k++)
@@ -571,8 +576,6 @@ 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.");
@@ -581,11 +584,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);
-
+ // 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++) {
@@ -610,11 +613,7 @@ 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("");
- }
+ PRINT_PART
/* If no derivative, double the smoothing length. */
if (wcount_dh == 0.0f) h_corr = p->h;
@@ -631,9 +630,11 @@ void runner_doghost(struct runner *r, struct cell *c) {
/* Did we get the right number density? */
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;
+
+ /* Ok, correct then */
+ p->h += h_corr;
+
+ /* And flag for another round of fun */
pid[redo] = pid[i];
redo += 1;
p->density.wcount = 0.0;
@@ -692,11 +693,9 @@ 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 */
@@ -704,8 +703,8 @@ void runner_doghost(struct runner *r, struct cell *c) {
count = redo;
if (count > 0) {
- //message("count=%d", count);
-
+ // message("count=%d", count);
+
/* Climb up the cell hierarchy. */
for (finger = c; finger != NULL; finger = finger->parent) {
@@ -784,6 +783,8 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
h = p->h;
ih = 1.0f / h;
+ PRINT_PART;
+
/* Drift... */
p->x[0] += xp->v_full[0] * dt;
p->x[1] += xp->v_full[1] * dt;
@@ -884,10 +885,6 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
m = p->mass;
x[0] = p->x[0], x[1] = p->x[1], x[2] = p->x[2];
- /* if (p->id == 0) */
- /* message("Kick ! t_beg=%f t_end=%f t_cur=%f", p->t_begin, p->t_end, */
- /* t_current); */
-
/* If particle needs to be kicked */
if (p->t_end <= t_current) {
@@ -926,6 +923,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
+
/* Kick particles in momentum space */
xp->v_full[0] += p->a[0] * dt;
xp->v_full[1] += p->a[1] * dt;