From 94dd38e63d26ddee40c99d891ab6c2819fdd0fb9 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sun, 28 Apr 2013 19:53:29 +0000
Subject: [PATCH] Added a new time step criterion. A particle can now not see
its volume (~h^3) change by more than a given constant factor over the course
of one time step.
Former-commit-id: c0a1cfcef2f3932b1c80ad7fc5c0001fa5711a1e
---
src/const.h | 1 +
src/runner.c | 11 +++++++----
2 files changed, 8 insertions(+), 4 deletions(-)
diff --git a/src/const.h b/src/const.h
index 98476d129f..249bcbf8f4 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 b6d82b5de4..f48133b500 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] );
--
GitLab