Skip to content
GitLab
Commit 628978c1 authored by Pedro Gonnet's avatar Pedro Gonnet
Browse files

several fixes while tracking down energy bugs. 


Former-commit-id: 9e7045749f44fad1d32ef1cafe9f73c7ae25d350
parent 5c13d956
Hide whitespace changes
Inline Side-by-side
src/cell.h
View file @ 628978c1
......@@ -20,6 +20,7 @@
/* Some constants. */
#define cell_sid_dt 13
/* The queue timers. */
enum {
cell_timer_none = 0,
......
src/const.h
View file @ 628978c1
......@@ -26,6 +26,7 @@
/* Time integration constants. */
#define const_cfl 0.3f
#define const_ln_max_h_change 0.231111721f /* Particle can't change volume by more than a factor of 2=1.26^3 over one time step */
#define const_max_u_change 0.1f
/* Neighbour search constants. */
#define const_eta_kernel 1.2348f /* Corresponds to 48 ngbs with the cubic spline kernel */
......
src/engine.c
View file @ 628978c1
......@@ -176,7 +176,7 @@ void engine_map_kick_first ( struct cell *c , void *data ) {
struct engine *e = (struct engine *)data;
float pdt, dt_step = e->dt_step, dt = e->dt, hdt = 0.5f*dt;
float dt_min, dt_max, h_max, dx, h2dx_max, dx_max;
float a[3], v[3], u, u_dt, h, h_dt, rho, v_old[3], w;
float a[3], v[3], u, u_dt, h, h_dt, v_old[3], w, rho;
double x[3], x_old[3];
struct part *restrict p;
struct xpart *restrict xp;
......@@ -236,14 +236,14 @@ void engine_map_kick_first ( struct cell *c , void *data ) {
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.05f )
/* if ( fabsf( w ) < 0.05f )
p->u = u *= 1.0f + w*( 1.0f + w*( -0.5f + w*1.0f/6.0f ) );
else
else */
p->u = u *= expf( w );
w = h_dt / h * dt;
if ( fabsf( w ) < 0.05f )
/* if ( fabsf( w ) < 0.05f )
p->h = h *= 1.0f + w*( 1.0f + w*( -0.5f + w*1.0f/6.0f ) );
else
else */
p->h = h *= expf( w );
h_max = fmaxf( h_max , h );
......@@ -259,11 +259,11 @@ void engine_map_kick_first ( struct cell *c , void *data ) {
/* Init fields for density calculation. */
if ( pdt > dt_step ) {
float w = -3.0f * h_dt / h * dt;
if ( fabsf( w ) < 0.1f )
/* if ( fabsf( w ) < 0.1f )
rho = p->rho *= 1.0f + w*( 1.0f + w*( -0.5f + w*(1.0f/6.0f - 1.0f/24.0*w ) ) );
else
else */
rho = p->rho *= expf( w );
p->force.POrho2 = u * ( const_gamma - 1.0f ) / ( rho + h * p->rho_dh / 3.0f );
p->force.POrho2 = u * ( const_gamma - 1.0f ) / ( rho * rho );
}
else {
p->density.wcount = 0.0f;
......@@ -407,7 +407,7 @@ void engine_single_density ( double *dim , long long int pid , struct part *__re
p.rho = ih * ih * ih * ( p.rho + p.mass*kernel_root );
p.rho_dh = p.rho_dh * ih * ih * ih * ih;
p.density.wcount = ( p.density.wcount + kernel_root ) * ( 4.0f / 3.0 * M_PI * kernel_gamma3 );
printf( "pairs_single: part %lli (h=%e) has wcount=%e, rho=%e, rho_dh=%e.\n" , p.id , p.h , p.density.wcount , p.rho , p.rho_dh );
printf( "pairs_single_density: part %lli (h=%e) has wcount=%e, rho=%e, rho_dh=%e.\n" , p.id , p.h , p.density.wcount , p.rho , p.rho_dh );
fflush(stdout);
}
......@@ -447,12 +447,17 @@ void engine_single_force ( double *dim , long long int pid , struct part *__rest
}
r2 = fdx[0]*fdx[0] + fdx[1]*fdx[1] + fdx[2]*fdx[2];
if ( r2 < p.h*p.h*kernel_gamma2 || r2 < parts[k].h*parts[k].h*kernel_gamma2 ) {
p.a[0] = 0.0f; p.a[1] = 0.0f; p.a[2] = 0.0f;
p.force.u_dt = 0.0f; p.force.h_dt = 0.0f; p.force.v_sig = 0.0f;
runner_iact_nonsym_force( r2 , fdx , p.h , parts[k].h , &p , &parts[k] );
double dvdr = ( (p.v[0]-parts[k].v[0])*fdx[0] + (p.v[1]-parts[k].v[1])*fdx[1] + (p.v[2]-parts[k].v[2])*fdx[2] ) / sqrt(r2);
printf( "pairs_single_force: part %lli and %lli interact (r=%.3e,dvdr=%.3e) with a=[%.3e,%.3e,%.3e], dudt=%.3e.\n" ,
p.id , parts[k].id , sqrt(r2) , dvdr , p.a[0] , p.a[1], p.a[2] , p.force.u_dt );
}
}
/* Dump the result. */
printf( "pairs_single: part %lli (h=%e) has a=[%.3e,%.3e,%.3e], udt=%e.\n" , p.id , p.h , p.a[0] , p.a[1] , p.a[2] , p.force.u_dt );
// printf( "pairs_single_force: part %lli (h=%e) has a=[%.3e,%.3e,%.3e], udt=%e.\n" , p.id , p.h , p.a[0] , p.a[1] , p.a[2] , p.force.u_dt );
fflush(stdout);
}
......@@ -478,18 +483,18 @@ void engine_step ( struct engine *e , int sort_queues ) {
TIMER_TIC2
/* Get the maximum dt. */
if ( e->policy & engine_policy_fixdt )
dt_step = FLT_MAX;
else {
if ( e->policy & engine_policy_multistep ) {
dt_step = 2.0f*dt;
for ( k = 0 ; k < 32 && (e->step & (1 << k)) == 0 ; k++ )
dt_step *= 2;
}
else
dt_step = FLT_MAX;
/* Set the maximum dt. */
e->dt_step = dt_step;
e->s->dt_step = dt_step;
printf( "engine_step: dt_step set to %.3e (dt=%.3e).\n" , dt_step , e->dt ); fflush(stdout);
// printf( "engine_step: dt_step set to %.3e (dt=%.3e).\n" , dt_step , e->dt ); fflush(stdout);
// printParticle( parts , 432626 );
......@@ -500,7 +505,7 @@ void engine_step ( struct engine *e , int sort_queues ) {
// for(k=0; k<10; ++k)
// printParticle(parts, k);
// printParticle( e->s->parts , 382557 , e->s->nr_parts );
// printParticle( e->s->parts , 3392063069037 , e->s->nr_parts );
/* Prepare the space. */
engine_prepare( e );
......@@ -514,7 +519,7 @@ void engine_step ( struct engine *e , int sort_queues ) {
}
}
// engine_single_density( e->s->dim , 7503966371841 , e->s->parts , e->s->nr_parts , e->s->periodic );
// engine_single_density( e->s->dim , 3392063069037 , e->s->parts , e->s->nr_parts , e->s->periodic );
/* Start the clock. */
TIMER_TIC_ND
......@@ -532,12 +537,13 @@ void engine_step ( struct engine *e , int sort_queues ) {
/* Stop the clock. */
TIMER_TOC(timer_runners);
// engine_single_force( e->s->dim , 7503966371841 , e->s->parts , e->s->nr_parts , e->s->periodic );
// engine_single_force( e->s->dim , 5497479945069 , e->s->parts , e->s->nr_parts , e->s->periodic );
// for(k=0; k<10; ++k)
// printParticle(parts, k);
// printParticle( parts , 432626 );
// printParticle( e->s->parts , 7503966371841 , e->s->nr_parts );
// printParticle( e->s->parts , 3392063069037 , e->s->nr_parts );
// printParticle( e->s->parts , 5497479945069 , e->s->nr_parts );
/* Collect the cell data from the second kick. */
for ( k = 0 ; k < e->s->nr_cells ; k++ ) {
......@@ -553,13 +559,17 @@ void engine_step ( struct engine *e , int sort_queues ) {
}
e->dt_min = dt_min;
e->dt_max = dt_max;
e->count_step = count;
e->ekin = ekin;
e->epot = epot;
// printParticle( e->s->parts , 382557 , e->s->nr_parts );
printf( "engine_step: dt_min/dt_max is %e/%e.\n" , dt_min , dt_max ); fflush(stdout);
printf( "engine_step: etot is %e (ekin=%e, epot=%e).\n" , ekin+epot , ekin , epot ); fflush(stdout);
printf( "engine_step: total momentum is [ %e , %e , %e ].\n" , mom[0] , mom[1] , mom[2] ); fflush(stdout);
printf( "engine_step: total angular momentum is [ %e , %e , %e ].\n" , ang[0] , ang[1] , ang[2] ); fflush(stdout);
// printf( "engine_step: dt_min/dt_max is %e/%e.\n" , dt_min , dt_max ); fflush(stdout);
// printf( "engine_step: etot is %e (ekin=%e, epot=%e).\n" , ekin+epot , ekin , epot ); fflush(stdout);
// printf( "engine_step: total momentum is [ %e , %e , %e ].\n" , mom[0] , mom[1] , mom[2] ); fflush(stdout);
// printf( "engine_step: total angular momentum is [ %e , %e , %e ].\n" , ang[0] , ang[1] , ang[2] ); fflush(stdout);
// printf( "engine_step: total entropic function is %e .\n", ent ); fflush(stdout);
printf( "engine_step: updated %i parts (dt_step=%.3e).\n" , count , dt_step ); fflush(stdout);
// printf( "engine_step: updated %i parts (dt_step=%.3e).\n" , count , dt_step ); fflush(stdout);
/* Increase the step. */
e->step += 1;
......@@ -576,18 +586,20 @@ void engine_step ( struct engine *e , int sort_queues ) {
e->dt *= 0.5;
while ( dt_min > 2*e->dt )
e->dt *= 2.0;
printf( "engine_step: dt_min=%.3e, adjusting time step to dt=%e.\n" , dt_min , e->dt );
// printf( "engine_step: dt_min=%.3e, adjusting time step to dt=%e.\n" , dt_min , e->dt );
}
else {
while ( dt_min < e->dt ) {
e->dt *= 0.5;
e->step *= 2;
printf( "engine_step: dt_min dropped below time step, adjusting to dt=%e.\n" , e->dt );
e->nullstep *= 2;
// printf( "engine_step: dt_min dropped below time step, adjusting to dt=%e.\n" , e->dt );
}
while ( dt_min > 2*e->dt && (e->step & 1) == 0 ) {
e->dt *= 2.0;
e->step /= 2;
printf( "engine_step: dt_min is larger than twice the time step, adjusting to dt=%e.\n" , e->dt );
e->nullstep /= 2;
// printf( "engine_step: dt_min is larger than twice the time step, adjusting to dt=%e.\n" , e->dt );
}
}
}
......
src/engine.h
View file @ 628978c1
......@@ -26,6 +26,7 @@
#define engine_policy_keep 4
#define engine_policy_block 8
#define engine_policy_fixdt 16
#define engine_policy_multistep 32
#define engine_queue_scale 1.2
......@@ -55,14 +56,20 @@ struct engine {
float dt_step;
/* The minimum dt over all particles in the system. */
float dt_min;
float dt_min, dt_max;
/* The system time step. */
float dt, dt_orig;
/* The system energies from the previous step. */
double ekin, epot;
/* The current step number. */
int step, nullstep;
/* The number of particles updated in the previous step. */
int count_step;
/* The current system time. */
float time;
......
src/runner.c
View file @ 628978c1
......@@ -519,7 +519,8 @@ void runner_dokick2 ( struct runner *r , struct cell *c ) {
float mom[3] = { 0.0f , 0.0f , 0.0f }, ang[3] = { 0.0f , 0.0f , 0.0f };
float x[3], v[3], u, h, pdt, m;
float dt_step = r->e->dt_step, dt = r->e->dt, hdt = 0.5f*dt;
float dt_cfl, dt_h_change;
float dt_cfl, dt_h_change, dt_u_change, dt_new;
float h_dt, u_dt;
struct part *p, *parts = c->parts;
struct xpart *xp;
......@@ -534,31 +535,36 @@ void runner_dokick2 ( struct runner *r , struct cell *c ) {
/* Get local copies of particle data. */
pdt = p->dt;
u_dt = p->force.u_dt;
h = p->h;
m = p->mass;
x[0] = p->x[0]; x[1] = p->x[1]; x[2] = p->x[2];
/* Scale the derivatives if they're freshly computed. */
if ( pdt <= dt_step ) {
p->force.h_dt *= h * 0.333333333f;
h_dt = p->force.h_dt *= h * 0.333333333f;
count += 1;
}
else
h_dt = p->force.h_dt;
/* Update the particle's time step. */
dt_cfl = const_cfl * h / p->force.v_sig;
dt_h_change = fabsf( const_ln_max_h_change * h / p->force.h_dt );
dt_h_change = ( h_dt != 0.0f ) ? fabsf( const_ln_max_h_change * h / h_dt ) : FLT_MAX;
dt_u_change = ( u_dt != 0.0f ) ? fabsf( const_max_u_change * p->u / u_dt ) : FLT_MAX;
dt_new = fminf( dt_cfl , fminf( dt_h_change , dt_u_change ) );
if ( pdt == 0.0f )
p->dt = pdt = fminf( dt_cfl , dt_h_change );
p->dt = pdt = dt_new;
else if ( pdt <= dt_step )
p->dt = pdt = fminf( fminf( dt_cfl , dt_h_change ) , 2.0f*pdt );
p->dt = pdt = fminf( dt_new , 2.0f*pdt );
else
p->dt = pdt = fminf( fminf( dt_cfl , dt_h_change ) , pdt );
p->dt = pdt = fminf( dt_new , pdt );
/* Update positions and energies at the half-step. */
p->v[0] = v[0] = xp->v_old[0] + hdt * p->a[0];
p->v[1] = v[1] = xp->v_old[1] + hdt * p->a[1];
p->v[2] = v[2] = xp->v_old[2] + hdt * p->a[2];
p->u = u = xp->u_old + hdt * p->force.u_dt;
p->u = u = xp->u_old + hdt * u_dt;
/* Get the smallest/largest dt. */
dt_min = fminf( dt_min , pdt );
......
src/runner_iact.h
View file @ 628978c1
......@@ -366,7 +366,7 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_force ( float r2
float wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, w, dvdr;
float mi, mj, POrho2i, POrho2j, rhoi, rhoj;
float v_sig, omega_ij, Pi_ij;
float dt_max;
// float dt_max;
float f;
int k;
......@@ -407,9 +407,9 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_force ( float r2
Pi_ij *= ( pi->force.balsara + pj->force.balsara );
/* Volker's modified viscosity */
dt_max = fmaxf(pi->dt, pj->dt);
if(dt_max > 0 && (wi_dr + wj_dr) < 0.)
Pi_ij = fminf( Pi_ij, 2.f * omega_ij / ( ( mi + mj ) * ( wi_dr + wj_dr ) * dt_max ) );
/* dt_max = fmaxf(pi->dt, pj->dt);
if( dt_max > 0 && (wi_dr + wj_dr) < 0. )
Pi_ij = fminf( Pi_ij, 2.f * omega_ij / ( ( mi + mj ) * ( wi_dr + wj_dr ) * dt_max ) ); */
/* Get the common factor out. */
w = ri * ( ( POrho2i * wi_dr + POrho2j * wj_dr ) + 0.25f * Pi_ij * ( wi_dr + wj_dr ) );
......@@ -600,14 +600,15 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_nonsym_force ( fl
float hi_inv, hi2_inv;
float hj_inv, hj2_inv;
float wi, wj, wi_dx, wj_dx, wi_dr, wj_dr, w, dvdr;
float mi, mj, POrho2i, POrho2j, rhoi, rhoj;
float /*mi,*/ mj, POrho2i, POrho2j, rhoi, rhoj;
float v_sig, omega_ij, Pi_ij;
float dt_max;
// float dt_max;
float f;
int k;
/* Get some values in local variables. */
mi = pi->mass; mj = pj->mass;
// mi = pi->mass;
mj = pj->mass;
rhoi = pi->rho; rhoj = pj->rho;
POrho2i = pi->force.POrho2;
POrho2j = pj->force.POrho2;
......@@ -643,9 +644,9 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_nonsym_force ( fl
Pi_ij *= ( pi->force.balsara + pj->force.balsara );
/* Volker's modified viscosity */
dt_max = fmaxf(pi->dt, pj->dt);
/* dt_max = fmaxf(pi->dt, pj->dt);
if(dt_max > 0 && (wi_dr + wj_dr) < 0.)
Pi_ij = fminf( Pi_ij, 2.f * omega_ij / ( ( mi + mj ) * ( wi_dr + wj_dr ) * dt_max ) );
Pi_ij = fminf( Pi_ij, 2.f * omega_ij / ( ( mi + mj ) * ( wi_dr + wj_dr ) * dt_max ) ); */
/* Get the common factor out. */
w = ri * ( ( POrho2i * wi_dr + POrho2j * wj_dr ) + 0.25f * Pi_ij * ( wi_dr + wj_dr ) );
......
src/space.c
View file @ 628978c1
......@@ -88,6 +88,7 @@ int space_marktasks ( struct space *s ) {
int k, nr_tasks = s->nr_tasks, *ind = s->tasks_ind;
struct task *t, *tasks = s->tasks;
float dt_step = s->dt_step;
struct cell *ci, *cj;
/* Run through the tasks and mark as skip or not. */
for ( k = 0 ; k < nr_tasks ; k++ ) {
......@@ -121,17 +122,22 @@ int space_marktasks ( struct space *s ) {
/* Set this task's skip. */
t->skip = ( t->ci->dt_min > dt_step && t->cj->dt_min > dt_step );
/* Local pointers. */
ci = t->ci;
cj = t->cj;
/* Too much particle movement? */
if ( t->tight &&
( t->ci->h2dx_max > t->ci->dmin || t->cj->h2dx_max > t->cj->dmin ) )
( ci->h2dx_max > ci->dmin || cj->h2dx_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 ) {
t->ci->sorts->flags |= (1 << t->flags);
t->ci->sorts->skip = 0;
t->cj->sorts->flags |= (1 << t->flags);
t->cj->sorts->skip = 0;
ci->sorts->flags |= (1 << t->flags);
ci->sorts->skip = 0;
cj->sorts->flags |= (1 << t->flags);
cj->sorts->skip = 0;
}
}
......
src/space.h
View file @ 628978c1
......@@ -29,6 +29,7 @@
#define space_dosub 1
#define space_stretch 1.05
#define space_maxtaskspercell 31
#define space_maxreldx 0.2f
/* Convert cell location to ID. */
......
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