Time-step calculation now includes a-factors

src/cosmology.c

@@ -137,6 +137,8 @@ void cosmology_update(struct cosmology *c, integertime_t ti_current) {
   c->a = a;
   c->a_inv = a_inv;
   c->a3_inv = a_inv * a_inv * a_inv;
+  c->a_factor_sig_vel =
+      pow(a, -1.5 * hydro_gamma_minus_one); /* a^{3*(1-gamma)/2} */
 
   /* Redshift */
   c->z = a_inv - 1.;
@@ -454,6 +456,9 @@ void cosmology_init_no_cosmo(const struct swift_params *params,
   c->H = 1.;
   c->a = 1.;
   c->a_inv = 1.;
+  c->a3_inv = 1.;
+  c->a_factor_sig_vel = 1.;
+
   c->a_dot = 0.;
   c->time = 0.;
   c->universe_age_at_present_day = 0.;

src/cosmology.h

@@ -41,6 +41,9 @@ struct cosmology {
   /*! Inverse cube of the current expansion factor of the Universe */
   double a3_inv;
 
+  /*! Power of the scale-factor used for signal-velocity conversion */
+  double a_factor_sig_vel;
+
   /*! Current redshift */
   double z;
 

src/gravity/Default/gravity.h

@@ -21,6 +21,7 @@
 #define SWIFT_DEFAULT_GRAVITY_H
 
 #include <float.h>
+#include "cosmology.h"
 #include "gravity_properties.h"
 #include "minmax.h"
 
@@ -67,7 +68,8 @@ __attribute__((always_inline)) INLINE static float gravity_get_potential(
  */
 __attribute__((always_inline)) INLINE static float
 gravity_compute_timestep_self(const struct gpart* const gp,
-                              const struct gravity_props* restrict grav_props) {
+                              const struct gravity_props* restrict grav_props,
+                              const struct cosmology* cosmo) {
 
   const float ac2 = gp->a_grav[0] * gp->a_grav[0] +
                     gp->a_grav[1] * gp->a_grav[1] +
@@ -78,7 +80,8 @@ gravity_compute_timestep_self(const struct gpart* const gp,
   const float epsilon = gravity_get_softening(gp);
 
   /* Note that 0.66666667 = 2. (from Gadget) / 3. (Plummer softening) */
-  const float dt = sqrtf(0.66666667f * grav_props->eta * epsilon * ac_inv);
+  const float dt =
+      sqrtf(0.66666667f * cosmo->a * grav_props->eta * epsilon * ac_inv);
 
   return dt;
 }

src/hydro/Gadget2/hydro.h

@@ -33,6 +33,7 @@
 
 #include "adiabatic_index.h"
 #include "approx_math.h"
+#include "cosmology.h"
 #include "dimension.h"
 #include "equation_of_state.h"
 #include "hydro_properties.h"
@@ -155,17 +156,19 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
  *
  * @param p Pointer to the particle data
  * @param xp Pointer to the extended particle data
+ * @param hydro_properties The constants used in the scheme
  *
  */
 __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
     const struct part *restrict p, const struct xpart *restrict xp,
-    const struct hydro_props *restrict hydro_properties) {
+    const struct hydro_props *restrict hydro_properties,
+    const struct cosmology *restrict cosmo) {
 
-  const float CFL_condition = hydro_properties->CFL_condition;
+  const float CFL = hydro_properties->CFL_condition;
 
   /* CFL condition */
-  const float dt_cfl =
-      2.f * kernel_gamma * CFL_condition * p->h / p->force.v_sig;
+  const float dt_cfl = 2.f * kernel_gamma * CFL * cosmo->a * p->h /
+                       (cosmo->a_factor_sig_vel * p->force.v_sig);
 
   return dt_cfl;
 }

src/timestep.h

@@ -37,10 +37,10 @@
  */
 __attribute__((always_inline)) INLINE static integertime_t
 make_integer_timestep(float new_dt, timebin_t old_bin, integertime_t ti_current,
-                      double timeBase_inv) {
+                      double time_base_inv) {
 
   /* Convert to integer time */
-  integertime_t new_dti = (integertime_t)(new_dt * timeBase_inv);
+  integertime_t new_dti = (integertime_t)(new_dt * time_base_inv);
 
   /* Current time-step */
   integertime_t current_dti = get_integer_timestep(old_bin);
@@ -51,7 +51,7 @@ make_integer_timestep(float new_dt, timebin_t old_bin, integertime_t ti_current,
 
   /* Put this timestep on the time line */
   integertime_t dti_timeline = max_nr_timesteps;
-  while (new_dti < dti_timeline) dti_timeline /= 2LL;
+  while (new_dti < dti_timeline) dti_timeline /= ((integertime_t)2);
   new_dti = dti_timeline;
 
   /* Make sure we are allowed to increase the timestep size */
@@ -70,16 +70,22 @@ make_integer_timestep(float new_dt, timebin_t old_bin, integertime_t ti_current,
 __attribute__((always_inline)) INLINE static integertime_t get_gpart_timestep(
     const struct gpart *restrict gp, const struct engine *restrict e) {
 
-  float new_dt = FLT_MAX;
+  float new_dt_self = FLT_MAX, new_dt_ext = FLT_MAX;
 
   if (e->policy & engine_policy_external_gravity)
-    new_dt =
-        min(new_dt, external_gravity_timestep(e->time, e->external_potential,
-                                              e->physical_constants, gp));
+    new_dt_ext = external_gravity_timestep(e->time, e->external_potential,
+                                           e->physical_constants, gp);
 
+  // MATTHIEU check accelerations a-factor!
   if (e->policy & engine_policy_self_gravity)
-    new_dt =
-        min(new_dt, gravity_compute_timestep_self(gp, e->gravity_properties));
+    new_dt_ext =
+        gravity_compute_timestep_self(gp, e->gravity_properties, e->cosmology);
+
+  /* Take the minimum of all */
+  float new_dt = min(new_dt_self, new_dt_ext);
+
+  /* Apply cosmology correction (H==1 if non-cosmological) */
+  new_dt *= e->cosmology->H;
 
   /* Limit timestep within the allowed range */
   new_dt = min(new_dt, e->dt_max);
@@ -89,7 +95,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_gpart_timestep(
 
   /* Convert to integer time */
   const integertime_t new_dti = make_integer_timestep(
-      new_dt, gp->time_bin, e->ti_current, e->timeBase_inv);
+      new_dt, gp->time_bin, e->ti_current, e->time_base_inv);
 
   return new_dti;
 }
@@ -106,7 +112,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
     const struct engine *restrict e) {
 
   /* Compute the next timestep (hydro condition) */
-  const float new_dt_hydro = hydro_compute_timestep(p, xp, e->hydro_properties);
+  const float new_dt_hydro =
+      hydro_compute_timestep(p, xp, e->hydro_properties, e->cosmology);
 
   /* Compute the next timestep (cooling condition) */
   float new_dt_cooling = FLT_MAX;
@@ -115,17 +122,19 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
                                       e->internal_units, p);
 
   /* Compute the next timestep (gravity condition) */
-  float new_dt_grav = FLT_MAX;
+  float new_dt_grav = FLT_MAX, new_dt_self_grav = FLT_MAX,
+        new_dt_ext_grav = FLT_MAX;
   if (p->gpart != NULL) {
 
     if (e->policy & engine_policy_external_gravity)
-      new_dt_grav = min(new_dt_grav, external_gravity_timestep(
-                                         e->time, e->external_potential,
-                                         e->physical_constants, p->gpart));
+      new_dt_ext_grav = external_gravity_timestep(
+          e->time, e->external_potential, e->physical_constants, p->gpart);
 
     if (e->policy & engine_policy_self_gravity)
-      new_dt_grav = min(new_dt_grav, gravity_compute_timestep_self(
-                                         p->gpart, e->gravity_properties));
+      new_dt_self_grav = gravity_compute_timestep_self(
+          p->gpart, e->gravity_properties, e->cosmology);
+
+    new_dt_grav = min(new_dt_self_grav, new_dt_ext_grav);
   }
 
   /* Final time-step is minimum of hydro and gravity */
@@ -139,6 +148,9 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
 
   new_dt = min(new_dt, dt_h_change);
 
+  /* Apply cosmology correction (H==1 if non-cosmological) */
+  new_dt *= e->cosmology->H;
+
   /* Limit timestep within the allowed range */
   new_dt = min(new_dt, e->dt_max);
   if (new_dt < e->dt_min)
@@ -147,7 +159,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
 
   /* Convert to integer time */
   const integertime_t new_dti = make_integer_timestep(
-      new_dt, p->time_bin, e->ti_current, e->timeBase_inv);
+      new_dt, p->time_bin, e->ti_current, e->time_base_inv);
 
   return new_dti;
 }
@@ -161,16 +173,25 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
 __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
     const struct spart *restrict sp, const struct engine *restrict e) {
 
-  float new_dt = star_compute_timestep(sp);
+  /* Stellar time-step */
+  float new_dt_star = star_compute_timestep(sp);
+
+  /* Gravity time-step */
+  float new_dt_self = FLT_MAX, new_dt_ext = FLT_MAX;
 
   if (e->policy & engine_policy_external_gravity)
-    new_dt = min(new_dt,
-                 external_gravity_timestep(e->time, e->external_potential,
-                                           e->physical_constants, sp->gpart));
+    new_dt_ext = external_gravity_timestep(e->time, e->external_potential,
+                                           e->physical_constants, sp->gpart);
 
   if (e->policy & engine_policy_self_gravity)
-    new_dt = min(new_dt, gravity_compute_timestep_self(sp->gpart,
-                                                       e->gravity_properties));
+    new_dt_self = gravity_compute_timestep_self(
+        sp->gpart, e->gravity_properties, e->cosmology);
+
+  /* Take the minimum of all */
+  float new_dt = min3(new_dt_star, new_dt_self, new_dt_ext);
+
+  /* Apply cosmology correction (H==1 if non-cosmological) */
+  new_dt *= e->cosmology->H;
 
   /* Limit timestep within the allowed range */
   new_dt = min(new_dt, e->dt_max);
@@ -180,7 +201,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
 
   /* Convert to integer time */
   const integertime_t new_dti = make_integer_timestep(
-      new_dt, sp->time_bin, e->ti_current, e->timeBase_inv);
+      new_dt, sp->time_bin, e->ti_current, e->time_base_inv);
 
   return new_dti;
 }