corrections to the time step selection.

Former-commit-id: ad42c06acc3d785ba5d38293a27e6f37b48b109c

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 */

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++ ) {

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];

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 ) );