From dd2a98704977b0c928f6d05aea9194b65a49d5cc Mon Sep 17 00:00:00 2001
From: Pedro Gonnet <pedro.gonnet@durham.ac.uk>
Date: Tue, 30 Apr 2013 15:20:07 +0000
Subject: [PATCH] corrections to the time step selection.
Former-commit-id: ad42c06acc3d785ba5d38293a27e6f37b48b109c
---
src/const.h | 2 +-
src/engine.c | 6 +++---
src/runner.c | 11 ++++++++---
src/runner_iact.h | 4 ++--
4 files changed, 14 insertions(+), 9 deletions(-)
diff --git a/src/const.h b/src/const.h
index 4530aa67d4..b686743e80 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 11a56b450b..ac906099f5 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 ec2fb6f5a8..bc9345c05a 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 faf3de134b..76f253048f 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 ) );
--
GitLab