diff --git a/src/const.h b/src/const.h index 98476d129fbfb9be7c1b48e041d38fa23010e417..249bcbf8f4a54d79142f8e45bdfbf1525a730ed3 100644 --- a/src/const.h +++ b/src/const.h @@ -25,6 +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 */ /* Neighbour search constants. */ #define const_nwneigh 48.f diff --git a/src/runner.c b/src/runner.c index b6d82b5de4dee5b35820561e3373bbe1533e88a7..f48133b500ea53d5d31197d0a037dfc1e4906837 100644 --- a/src/runner.c +++ b/src/runner.c @@ -519,6 +519,7 @@ 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; struct part *p, *parts = c->parts; struct xpart *xp; @@ -544,7 +545,9 @@ void runner_dokick2 ( struct runner *r , struct cell *c ) { } /* Update the particle's time step. */ - p->dt = pdt = const_cfl * h / p->force.v_sig; + 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); /* Update positions and energies at the half-step. */ p->v[0] = v[0] = xp->v_old[0] + hdt * p->a[0]; @@ -561,9 +564,9 @@ void runner_dokick2 ( struct runner *r , struct cell *c ) { epot += m * u; /* Collect momentum */ - mom[0] += m * p->v[0]; - mom[1] += m * p->v[1]; - mom[2] += m * p->v[2]; + mom[0] += m * v[0]; + mom[1] += m * v[1]; + mom[2] += m * v[2]; /* Collect angular momentum */ ang[0] += m * ( x[1]*v[2] - x[2]*v[1] );