Commit 0a89ef0f authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Can now do the initialisation and one time step. Still freezes in the second time step.

parent 5446b651
......@@ -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);
}
......
......@@ -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;
......
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