diff --git a/src/const.h b/src/const.h index 4530aa67d4295eaa45dd5fd256ddebb9f4bc8e3a..b686743e806530ee04ff5528b2cf2d0fa32ae4ca 100644 --- a/src/const.h +++ b/src/const.h @@ -25,7 +25,7 @@ /* Time integration constants. */ #define const_cfl 0.3f -#define const_ln_max_h_change log(1.26f) /* Particle can't change volume by more than a factor of 2=1.26^3 over one time step */ +#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 */ /* Neighbour search constants. */ #define const_eta_kernel 1.1272f /* Corresponds to 48 ngbs with the cubic spline kernel */ diff --git a/src/engine.c b/src/engine.c index 11a56b450b85826474d0123895e3eaccd1179234..ac906099f53857e1dbff7597b11dd476f5ddc9b7 100644 --- a/src/engine.c +++ b/src/engine.c @@ -514,7 +514,7 @@ void engine_step ( struct engine *e , int sort_queues ) { } } - // engine_single_density( e->s->dim , 494748 , e->s->parts , e->s->nr_parts , e->s->periodic ); + // engine_single_density( e->s->dim , 7503966371841 , e->s->parts , e->s->nr_parts , e->s->periodic ); /* Start the clock. */ TIMER_TIC_ND @@ -532,12 +532,12 @@ void engine_step ( struct engine *e , int sort_queues ) { /* Stop the clock. */ TIMER_TOC(timer_runners); - // engine_single_force( e->s->dim , 494748 , e->s->parts , e->s->nr_parts , e->s->periodic ); + // engine_single_force( e->s->dim , 7503966371841 , 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 , 494748 , e->s->nr_parts ); + // printParticle( e->s->parts , 7503966371841 , e->s->nr_parts ); /* Collect the cell data from the second kick. */ for ( k = 0 ; k < e->s->nr_cells ; k++ ) { diff --git a/src/runner.c b/src/runner.c index ec2fb6f5a87547f36d957635ddcd70e61540f3dc..bc9345c05a79f99f27ecd9bb2ea92c4c8feffb6f 100644 --- a/src/runner.c +++ b/src/runner.c @@ -545,9 +545,14 @@ void runner_dokick2 ( struct runner *r , struct cell *c ) { } /* Update the particle's time step. */ - dt_cfl = const_cfl * h / p->force.v_sig; - dt_h_change = const_ln_max_h_change * h / p->force.h_dt; - p->dt = pdt = fminf(dt_cfl, dt_h_change); + dt_cfl = const_cfl * h / p->force.v_sig; + dt_h_change = fabsf( const_ln_max_h_change * h / p->force.h_dt ); + if ( pdt == 0.0f ) + p->dt = pdt = fminf( dt_cfl , dt_h_change ); + else if ( pdt <= dt_step ) + p->dt = pdt = fminf( fminf( dt_cfl , dt_h_change ) , 2.0f*pdt ); + else + p->dt = pdt = fminf( fminf( dt_cfl , dt_h_change ) , pdt ); /* Update positions and energies at the half-step. */ p->v[0] = v[0] = xp->v_old[0] + hdt * p->a[0]; diff --git a/src/runner_iact.h b/src/runner_iact.h index faf3de134b4194d339a2f7fc6624a36ceec14587..76f253048f4387df1eb1a3b70c391e404ca2ac0c 100644 --- a/src/runner_iact.h +++ b/src/runner_iact.h @@ -409,7 +409,7 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_force ( float r2 /* 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 ) ); + 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 ) ); @@ -645,7 +645,7 @@ __attribute__ ((always_inline)) INLINE static void runner_iact_nonsym_force ( fl /* 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 ) ); + 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 ) );