diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 8fcfe6490c539732b885a9c102f011aebb0aaeff..0906c2e344671bb58cc7fbdb0a5a86028e083578 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -28,9 +28,12 @@
__attribute__((always_inline))
INLINE static float gravity_compute_timestep(struct gpart* gp) {
+ float dt = FLT_MAX;
- /* Currently no limit is imposed */
- return FLT_MAX;
+#ifdef EXTERNAL_POTENTIAL_POINTMASS
+ dt = fminf(dt, external_gravity_pointmass_timestep(gp));
+#endif
+ return dt;
}
/**
@@ -66,6 +69,7 @@ __attribute__((always_inline)) INLINE static void external_gravity(struct gpart
#ifdef EXTERNAL_POTENTIAL_POINTMASS
external_gravity_pointmass(g);
#endif
+}
/**
* @brief Kick the additional variables
@@ -76,3 +80,6 @@ __attribute__((always_inline)) INLINE static void external_gravity(struct gpart
*/
__attribute__((always_inline)) INLINE static void gravity_kick_extra(
struct gpart* gp, float dt, float half_dt) {}
+
+__attribute__((always_inline)) INLINE static void gravity_end_force(
+ struct gpart* gp) {}
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 47aadd19585a71d04443403245cc76e28f46d80b..6a33f786cd2ee1bb6fdbba20e544fa4b9d2dd556 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -31,9 +31,6 @@ struct gpart {
/* Particle acceleration. */
float a_grav[3];
- /* Gravity acceleration */
- float a_grav_external[3];
-
/* Particle mass. */
float mass;
diff --git a/src/gravity/Default/potentials.h b/src/gravity/Default/potentials.h
index 39cde2d1dd50631615ab5474b15d85b4302b0746..9930b66ec246a18cdb9b71414f919baeaeb274c4 100644
--- a/src/gravity/Default/potentials.h
+++ b/src/gravity/Default/potentials.h
@@ -25,7 +25,7 @@ __attribute__((always_inline)) INLINE static float external_gravity_pointmass_ti
const float dota_y = const_G * External_Potential_Mass * rinv * rinv * rinv * (-g->v_full[1] + 3.f * rinv * rinv * drdv * dy);
const float dota_z = const_G * External_Potential_Mass * rinv * rinv * rinv * (-g->v_full[2] + 3.f * rinv * rinv * drdv * dz);
const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
- const float a_2 = g->a_grav_external[0] * g->a_grav_external[0] + g->a_grav_external[1] * g->a_grav_external[1] + g->a_grav_external[2] * g->a_grav_external[2];
+ const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] + g->a_grav[2] * g->a_grav[2];
return 0.03f * sqrtf(a_2/dota_2);
}
@@ -38,9 +38,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_pointmass(str
const float rinv = 1.f / sqrtf(dx*dx + dy*dy + dz*dz);
- g->a_grav_external[0] += - const_G * External_Potential_Mass * dx * rinv * rinv * rinv;
- g->a_grav_external[1] += - const_G * External_Potential_Mass * dy * rinv * rinv * rinv;
- g->a_grav_external[2] += - const_G * External_Potential_Mass * dz * rinv * rinv * rinv;
+ g->a_grav[0] += - const_G * External_Potential_Mass * dx * rinv * rinv * rinv;
+ g->a_grav[1] += - const_G * External_Potential_Mass * dy * rinv * rinv * rinv;
+ g->a_grav[2] += - const_G * External_Potential_Mass * dz * rinv * rinv * rinv;
}
diff --git a/src/runner.c b/src/runner.c
index 3341f3213261ff8d81625e9dc2103148b169bd32..b6da14bfe34b1ec3415fcdaabec227d130e995d6 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -139,16 +139,16 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
const float rinv = 1.f / sqrtf(dx*dx + dy*dy + dz*dz);
/* /\* current acceleration *\/ */
- /* g->a_grav_external[0] = - const_G * External_Potential_Mass * dx * rinv * rinv * rinv; */
- /* g->a_grav_external[1] = - const_G * External_Potential_Mass * dy * rinv * rinv * rinv; */
- /* g->a_grav_external[2] = - const_G * External_Potential_Mass * dz * rinv * rinv * rinv; */
+ /* g->a_grav[0] = - const_G * External_Potential_Mass * dx * rinv * rinv * rinv; */
+ /* g->a_grav[1] = - const_G * External_Potential_Mass * dy * rinv * rinv * rinv; */
+ /* g->a_grav[2] = - const_G * External_Potential_Mass * dz * rinv * rinv * rinv; */
const int current_dti = g->ti_end - g->ti_begin;
const float dt = 0.5f * current_dti * r->e->timeBase;
- const float vx = g->v_full[0] + dt * g->a_grav_external[0];
- const float vy = g->v_full[1] + dt * g->a_grav_external[1];
- const float vz = g->v_full[2] + dt * g->a_grav_external[2];
+ const float vx = g->v_full[0] + dt * g->a_grav[0];
+ const float vy = g->v_full[1] + dt * g->a_grav[1];
+ const float vz = g->v_full[2] + dt * g->a_grav[2];
/* E/L */
E = 0.5 * ((vx*vx) + (vy*vy) + (vz*vz)) - const_G * External_Potential_Mass * rinv;
@@ -158,7 +158,7 @@ void runner_dograv_external(struct runner *r, struct cell *c) {
if(abs(g->id) == 1) {
float v2 = vx*vx + vy*vy + vz*vz;
float fg = const_G * External_Potential_Mass * rinv;
- float fga = sqrtf((g->a_grav_external[0]*g->a_grav_external[0]) + (g->a_grav_external[1]*g->a_grav_external[1]) + (g->a_grav_external[2]*g->a_grav_external[2])) * dr;
+ float fga = sqrtf((g->a_grav[0]*g->a_grav[0]) + (g->a_grav[1]*g->a_grav[1]) + (g->a_grav[2]*g->a_grav[2])) * dr;
//message("grav_external time= %f\t V_c^2= %f GM/r= %f E= %f L[2]= %f x= %f y= %f vx= %f vy= %f\n", r->e->time, v2, fg, E, L[2], g->x[0], g->x[1], vx, vy);
message("%f\t %f %f %f %f %f %f %f %f %f %f %f %f\n", r->e->time, g->tx, g->tv, dt, v2, fg, fga, E, L[2], g->x[0], g->x[1], vx, vy);
// message(" G=%e M=%e\n", const_G, External_Potential_Mass);
@@ -955,6 +955,10 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
const double dt = (ti_end - ti_start) * timeBase;
const double half_dt = (ti_end - gp->ti_end) * timeBase;
+ /* Move particle forward in time */
+ gp->ti_begin = gp->ti_end;
+ gp->ti_end = gp->ti_begin + new_dti;
+
/* Kick particles in momentum space */
gp->v_full[0] += gp->a_grav[0] * dt;
gp->v_full[1] += gp->a_grav[1] * dt;
@@ -966,6 +970,11 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* Number of updated g-particles */
g_updated++;
}
+
+ /* Minimal time for next end of time-step */
+ ti_end_min = min(gp->ti_end, ti_end_min);
+ ti_end_max = max(gp->ti_end, ti_end_max);
+
}
/* Now do the hydro ones... */