Commit e0bee388 authored by Pedro Gonnet's avatar Pedro Gonnet
Browse files

fix multiple time steps, i.e. mindt should now work, multistep is still not quite right.


Former-commit-id: fcde932ed0b4842501793f8d7d1c71a50522554a
parent 565675f2
......@@ -46,10 +46,8 @@
cell distance over the cell width. */
// #define const_G 6.67384e-8f /* Gravitational constant. */
#define const_G 6.672e-8f /* Gravitational constant. */
// #define const_epsilon 0.0014f /* Gravity blending distance. */
// #define const_iepsilon 714.285714286f /* Inverse gravity blending distance. */
#define const_epsilon 1e-20
#define const_iepsilon 1e20
#define const_epsilon 0.0014f /* Gravity blending distance. */
#define const_iepsilon 714.285714286f /* Inverse gravity blending distance. */
#define const_iepsilon2 (const_iepsilon*const_iepsilon)
#define const_iepsilon3 (const_iepsilon2*const_iepsilon)
#define const_iepsilon4 (const_iepsilon2*const_iepsilon2)
......
......@@ -1063,10 +1063,6 @@ void engine_maketasks ( struct engine *e ) {
super cells. */
for ( k = 0 ; k < nr_cells ; k++ )
engine_mkghosts( e , &cells[k] , NULL );
/* if ( e->policy & engine_policy_fixdt )
space_map_cells_pre( s , 1 , &scheduler_map_mkghosts_nokick1 , sched );
else
space_map_cells_pre( s , 1 , &scheduler_map_mkghosts , sched ); */
/* Run through the tasks and make force tasks for each density task.
Each force task depends on the cell ghosts and unlocks the kick2 task
......@@ -1395,14 +1391,14 @@ void engine_prepare ( struct engine *e ) {
if ( rebuild ) {
// tic = getticks();
engine_rebuild( e );
// message( "engine_rebuild took %.3f ms." , (double)(getticks() - tic)/CPU_TPS*1000 );
// message( "engine_rebuild took %.3f ms." , (double)(getticks() - tic)/CPU_TPS*1000 );
}
/* Re-rank the tasks every now and then. */
if ( e->tasks_age % engine_tasksreweight == 1 ) {
// tic = getticks();
scheduler_reweight( &e->sched );
// message( "scheduler_reweight took %.3f ms." , (double)(getticks() - tic)/CPU_TPS*1000 );
// message( "scheduler_reweight took %.3f ms." , (double)(getticks() - tic)/CPU_TPS*1000 );
}
e->tasks_age += 1;
......@@ -2095,7 +2091,6 @@ void engine_init ( struct engine *e , struct space *s , float dt , int nr_thread
scheduler_reset( &e->sched , s->tot_cells );
for ( k = 0 ; k < s->nr_cells ; k++ )
s->cells[k].kick1 = scheduler_addtask( &e->sched , task_type_kick1 , task_subtype_none , 0 , 0 , &s->cells[k] , NULL , 0 );
// space_map_cells_pre( e->s , 1 , &scheduler_map_mkkick1 , &e->sched );
scheduler_ranktasks( &e->sched );
/* Allocate and init the threads. */
......
......@@ -748,11 +748,11 @@ void runner_doghost ( struct runner *r , struct cell *c ) {
void runner_dokick2 ( struct runner *r , struct cell *c ) {
int j, k, count = 0, nr_parts = c->count;
float dt_curr, hdt_curr, dt_min = FLT_MAX, dt_max = 0.0f;
float dt_min = FLT_MAX, dt_max = 0.0f;
double ekin = 0.0, epot = 0.0;
float mom[3] = { 0.0f , 0.0f , 0.0f }, ang[3] = { 0.0f , 0.0f , 0.0f };
float x[3], v_hdt[3], u_hdt, h, pdt, m;
float dt_step = r->e->dt_step, dt = r->e->dt, idt;
float dt_step = r->e->dt_step, dt = r->e->dt, hdt, idt;
float dt_cfl, dt_h_change, dt_u_change, dt_new;
float h_dt, u_dt;
struct part *restrict p, *restrict parts = c->parts;
......@@ -762,6 +762,7 @@ void runner_dokick2 ( struct runner *r , struct cell *c ) {
/* Init idt to avoid compiler stupidity. */
idt = ( dt > 0 ) ? 1.0f / dt : 0.0f;
hdt = dt / 2;
/* Loop over the particles and kick them. */
__builtin_prefetch( &parts[0] , 0 , 1 );
......@@ -811,23 +812,15 @@ void runner_dokick2 ( struct runner *r , struct cell *c ) {
else
p->dt = pdt = fminf( dt_new , 2.0f*pdt );
/* Get the particle-specific time step. */
dt_curr = xp->dt_curr;
hdt_curr = 0.5f * dt_curr;
/* Update positions and energies at the full step. */
p->v[0] = v_hdt[0] + hdt_curr * p->a[0];
p->v[1] = v_hdt[1] + hdt_curr * p->a[1];
p->v[2] = v_hdt[2] + hdt_curr * p->a[2];
p->u = u_hdt + hdt_curr * u_dt;
xp->v_hdt[0] = ( v_hdt[0] += dt_curr * p->a[0] );
xp->v_hdt[1] = ( v_hdt[1] += dt_curr * p->a[1] );
xp->v_hdt[2] = ( v_hdt[2] += dt_curr * p->a[2] );
xp->u_hdt = ( u_hdt += dt_curr * u_dt );
p->v[0] = v_hdt[0] + hdt * p->a[0];
p->v[1] = v_hdt[1] + hdt * p->a[1];
p->v[2] = v_hdt[2] + hdt * p->a[2];
p->u = u_hdt + hdt * u_dt;
/* Set the new particle-specific time step. */
if ( dt > 0.0f ) {
dt_curr = dt;
float dt_curr = dt;
j = (int)( pdt * idt );
while ( j > 1 ) {
dt_curr *= 2.0f;
......@@ -890,7 +883,7 @@ void runner_dokick1 ( struct runner *r , struct cell *c ) {
int j, k;
struct engine *e = r->e;
float pdt, dt_step = e->dt_step, dt = e->dt;
float pdt, dt_step = e->dt_step, dt = e->dt, hdt = dt/2;
float dt_min, dt_max, h_max, dx, dx_max;
float a[3], v[3], u, u_dt, h, h_dt, w, rho;
double x[3], x_old[3];
......@@ -940,6 +933,12 @@ void runner_dokick1 ( struct runner *r , struct cell *c ) {
dt_min = fminf( dt_min , pdt );
dt_max = fmaxf( dt_max , pdt );
/* Update the half-step velocities from the current velocities. */
xp->v_hdt[0] = v[0] + hdt * a[0];
xp->v_hdt[1] = v[1] + hdt * a[1];
xp->v_hdt[2] = v[2] + hdt * a[2];
xp->u_hdt = u + hdt * u_dt;
/* Move the particles with the velocities at the half-step. */
p->x[0] = x[0] += dt * xp->v_hdt[0];
p->x[1] = x[1] += dt * xp->v_hdt[1];
......@@ -950,9 +949,9 @@ void runner_dokick1 ( struct runner *r , struct cell *c ) {
dx_max = fmaxf( dx_max , dx );
/* Update positions and energies at the half-step. */
p->v[0] = v[0] += dt * a[0];
p->v[1] = v[1] += dt * a[1];
p->v[2] = v[2] += dt * a[2];
p->v[0] = v[0] + dt * a[0];
p->v[1] = v[1] + dt * a[1];
p->v[2] = v[2] + dt * a[2];
w = u_dt / u * dt;
if ( fabsf( w ) < 0.01f )
p->u = u *= 1.0f + w*( 1.0f + w*( 0.5f + w*( 1.0f/6.0f + 1.0f/24.0f*w ) ) );
......
......@@ -99,29 +99,6 @@ void scheduler_addunlock ( struct scheduler *s , struct task *ta , struct task *
}
/**
* @brief Mapping function to append a ghost task to each cell.
*
* A kick1-task is appended to each cell.
*/
void scheduler_map_mkkick1 ( struct cell *c , void *data ) {
struct scheduler *s = (struct scheduler *)data;
struct cell *finger;
/* Append a kick1 task and make sure the parent depends on it. */
c->kick1 = scheduler_addtask( s , task_type_kick1 , task_subtype_none , 0 , 0 , c , NULL , 0 );
if ( c->parent != NULL )
scheduler_addunlock( s , c->kick1 , c->parent->kick1 );
/* Set a bogus super cell. */
for ( finger = c ; finger->parent != NULL ; finger = finger->parent );
c->super = finger;
}
/**
* @brief Split tasks that may be too large.
*
......@@ -837,7 +814,7 @@ void scheduler_start ( struct scheduler *s , unsigned int mask ) {
for ( k = nr_tasks - 1 ; k >= 0 ; k-- ) {
t = &tasks[ tid[k] ];
t->wait = 0;
t->rid = -1;
t->rid = -1;
if ( !( (1 << t->type) & mask ) || t->skip )
continue;
for ( j = 0 ; j < t->nr_unlock_tasks ; j++ )
......@@ -845,17 +822,20 @@ void scheduler_start ( struct scheduler *s , unsigned int mask ) {
}
// message( "waiting tasks took %.3f ms." , (double)( getticks() - tic ) / CPU_TPS * 1000 );
/* Don't enqueue link tasks directly. */
mask &= ~(1 << task_type_link);
/* Loop over the tasks and enqueue whoever is ready. */
// tic = getticks();
for ( k = 0 ; k < nr_tasks ; k++) {
t = &tasks[ tid[k] ];
if ( ( (1 << t->type) & mask ) && !t->skip ) {
if ( t->wait == 0 ) {
scheduler_enqueue( s , t );
pthread_cond_broadcast( &s->sleep_cond );
}
else
break;
scheduler_enqueue( s , t );
pthread_cond_broadcast( &s->sleep_cond );
}
else
break;
}
}
// message( "enqueueing tasks took %.3f ms." , (double)( getticks() - tic ) / CPU_TPS * 1000 );
......
......@@ -81,7 +81,6 @@ void scheduler_ranktasks ( struct scheduler *s );
void scheduler_reweight ( struct scheduler *s );
struct task *scheduler_addtask ( struct scheduler *s , int type , int subtype , int flags , int wait , struct cell *ci , struct cell *cj , int tight );
void scheduler_splittasks ( struct scheduler *s );
void scheduler_map_mkkick1 ( struct cell *c , void *data );
struct task *scheduler_done ( struct scheduler *s , struct task *t );
struct task *scheduler_unlock ( struct scheduler *s , struct task *t );
void scheduler_addunlock ( struct scheduler *s , struct task *ta , struct task *tb );
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