Replace ti_begin and ti_end by a time bin. Still need to fix cooling(),...

Replace ti_begin and ti_end by a time bin. Still need to fix cooling(), do_recv() statistics() and propagate types everywhere.

examples/SedovBlast_3D/makeIC.py

@@ -54,6 +54,7 @@ u[:] = P0 / (rho0 * (gamma - 1))
 # Make the central particles detonate
 index = argsort(r)
 u[index[0:N_inject]] = E0 / (N_inject * m[0])
+print ids[index[0:N_inject]]
 
 #--------------------------------------------------
 

src/active.h

@@ -101,14 +101,18 @@ __attribute__((always_inline)) INLINE static int cell_is_all_active(
 __attribute__((always_inline)) INLINE static int part_is_active(
     const struct part *p, const struct engine *e) {
 
+  const integertime_t ti_current = e->ti_current;
+  const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
 #ifdef SWIFT_DEBUG_CHECKS
-  if (p->ti_end < e->ti_current)
-    error("particle in an impossible time-zone! p->ti_end=%d e->ti_current=%d",
-          p->ti_end, e->ti_current);
+  if (ti_end < ti_current)
+    error(
+        "particle in an impossible time-zone! p->ti_end=%lld "
+        "e->ti_current=%lld",
+        ti_end, ti_current);
 #endif
 
-  // return (p->ti_end == e->ti_current);
-  return (p->time_bin < get_max_active_bin(e->ti_current));
+  return (ti_end == ti_current);
 }
 
 /**
@@ -120,15 +124,18 @@ __attribute__((always_inline)) INLINE static int part_is_active(
 __attribute__((always_inline)) INLINE static int gpart_is_active(
     const struct gpart *gp, const struct engine *e) {
 
+  const integertime_t ti_current = e->ti_current;
+  const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
+
 #ifdef SWIFT_DEBUG_CHECKS
-  if (gp->ti_end < e->ti_current)
+  if (ti_end < ti_current)
     error(
-        "g-particle in an impossible time-zone! gp->ti_end=%d e->ti_current=%d",
-        gp->ti_end, e->ti_current);
+        "g-particle in an impossible time-zone! gp->ti_end=%lld "
+        "e->ti_current=%lld",
+        ti_end, ti_current);
 #endif
 
-  // return (gp->ti_end == e->ti_current);
-  return (gp->time_bin < get_max_active_bin(e->ti_current));
+  return (ti_end == ti_current);
 }
 
 #endif /* SWIFT_ACTIVE_H */

src/engine.c

@@ -2643,6 +2643,10 @@ void engine_step(struct engine *e) {
   engine_launch(e, e->nr_threads);
   TIMER_TOC(timer_runners);
 
+  for (size_t i = 0; i < e->s->nr_parts; ++i) {
+    if (e->s->parts[i].time_bin == 0) error("Particle in bin 0");
+  }
+
   /* Save some statistics */
   if (e->time - e->timeLastStatistics >= e->deltaTimeStatistics) {
     engine_print_stats(e);

src/gravity/Default/gravity.h

@@ -53,7 +53,7 @@ gravity_compute_timestep_self(const struct gpart* const gp) {
 __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
     struct gpart* gp) {
 
-  gp->time_bin = -1;
+  gp->time_bin = 0;
   gp->epsilon = 0.;  // MATTHIEU
 }
 

src/gravity/Default/gravity_debug.h

@@ -24,10 +24,9 @@ __attribute__((always_inline)) INLINE static void gravity_debug_particle(
   printf(
       "x=[%.3e,%.3e,%.3e], "
       "v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
-      "mass=%.3e t_begin=%d, t_end=%d\n",
+      "mass=%.3e time_bin=%d\n",
       p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
-      p->a_grav[0], p->a_grav[1], p->a_grav[2], p->mass, p->ti_begin,
-      p->ti_end);
+      p->a_grav[0], p->a_grav[1], p->a_grav[2], p->mass, p->time_bin);
 }
 
 #endif /* SWIFT_DEFAULT_GRAVITY_DEBUG_H */

src/hydro/Gadget2/hydro.h

@@ -212,7 +212,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
 __attribute__((always_inline)) INLINE static void hydro_first_init_part(
     struct part *restrict p, struct xpart *restrict xp) {
 
-  p->time_bin = -1;
+  p->time_bin = 0;
   xp->v_full[0] = p->v[0];
   xp->v_full[1] = p->v[1];
   xp->v_full[2] = p->v[2];
@@ -385,8 +385,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     p->rho *= expf(w2);
 
   /* Drift the pressure */
-  const integertime_t ti_begin = get_integer_time_begin(t1, p->time_bin);
-  const integertime_t ti_end = get_integer_time_end(t1, p->time_bin);
+  const integertime_t ti_begin = get_integer_time_begin(t0, p->time_bin);
+  const integertime_t ti_end = get_integer_time_end(t0, p->time_bin);
   const float dt_entr = (t1 - (ti_begin + ti_end) / 2) * timeBase;
   const float pressure = hydro_get_pressure(p, dt_entr);
 

src/hydro/Gadget2/hydro_debug.h

@@ -27,14 +27,14 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       "h=%.3e, wcount=%.3f, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, "
       "P=%.3e, P_over_rho2=%.3e, S=%.3e, dS/dt=%.3e, c=%.3e\n"
       "divV=%.3e, rotV=[%.3e,%.3e,%.3e], balsara=%.3e \n "
-      "v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
+      "v_sig=%e dh/dt=%.3e time_bin=%d\n",
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
       xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
       p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh,
       p->rho, hydro_get_pressure(p, 0.), p->force.P_over_rho2, p->entropy,
       p->entropy_dt, p->force.soundspeed, p->density.div_v, p->density.rot_v[0],
       p->density.rot_v[1], p->density.rot_v[2], p->force.balsara,
-      p->force.v_sig, p->force.h_dt, p->ti_begin, p->ti_end);
+      p->force.v_sig, p->force.h_dt, p->time_bin);
 }
 
 #endif /* SWIFT_GADGET2_HYDRO_DEBUG_H */

src/kick.h

@@ -32,20 +32,27 @@
  *
  * @param gp The #gpart to kick.
  * @param new_dti The (integer) time-step for this kick.
+ * @param ti_current The current (integer) time.
  * @param timeBase The minimal allowed time-step size.
  */
 __attribute__((always_inline)) INLINE static void kick_gpart(
-    struct gpart *restrict gp, int new_dti, double timeBase) {
+    struct gpart *restrict gp, integertime_t new_dti, integertime_t ti_current,
+    double timeBase) {
 
   /* Compute the time step for this kick */
-  const int ti_start = (gp->ti_begin + gp->ti_end) / 2;
-  const int ti_end = gp->ti_end + new_dti / 2;
+  const integertime_t old_ti_begin =
+      get_integer_time_begin(ti_current, gp->time_bin);
+  const integertime_t old_ti_end =
+      get_integer_time_end(ti_current, gp->time_bin);
+  const int ti_start = (old_ti_begin + old_ti_end) / 2;
+  const int ti_end = old_ti_end + new_dti / 2;
   const float dt = (ti_end - ti_start) * timeBase;
-  const float half_dt = (ti_end - gp->ti_end) * timeBase;
+  const float half_dt = (ti_end - old_ti_end) * timeBase;
 
   /* Move particle forward in time */
-  gp->ti_begin = gp->ti_end;
-  gp->ti_end = gp->ti_begin + new_dti;
+  // gp->ti_begin = gp->ti_end;
+  // gp->ti_end = gp->ti_begin + new_dti;
+  gp->time_bin = get_time_bin(new_dti);
 
   /* Kick particles in momentum space */
   gp->v_full[0] += gp->a_grav[0] * dt;
@@ -62,26 +69,33 @@ __attribute__((always_inline)) INLINE static void kick_gpart(
  * @param p The #part to kick.
  * @param xp The #xpart of the particle.
  * @param new_dti The (integer) time-step for this kick.
+ * @param ti_current The current (integer) time.
  * @param timeBase The minimal allowed time-step size.
  */
 __attribute__((always_inline)) INLINE static void kick_part(
     struct part *restrict p, struct xpart *restrict xp, int new_dti,
-    double timeBase) {
+    integertime_t ti_current, double timeBase) {
 
   /* Compute the time step for this kick */
-  const integertime_t ti_begin = get_integer_time_begin(t1, p->time_bin);
-  const integertime_t ti_end = get_integer_time_end(t1, p->time_bin);
-  const int ti_start = (ti_begin + ti_end) / 2;
-  const int ti_end = ti_end + new_dti / 2;
+  const integertime_t old_ti_begin =
+      get_integer_time_begin(ti_current, p->time_bin);
+  const integertime_t old_ti_end =
+      get_integer_time_end(ti_current, p->time_bin);
+  const int ti_start = (old_ti_begin + old_ti_end) / 2;
+  const int ti_end = old_ti_end + new_dti / 2;
   const float dt = (ti_end - ti_start) * timeBase;
-  const float half_dt = (ti_end - p->ti_end) * timeBase;
+  const float half_dt = (ti_end - old_ti_end) * timeBase;
 
   /* Move particle forward in time */
-  p->ti_begin = p->ti_end;
-  p->ti_end = p->ti_begin + new_dti;
+  // p->ti_begin = p->ti_end;
+  // p->ti_end = p->ti_begin + new_dti;
+  p->time_bin = get_time_bin(new_dti);
+  // if (p->id == 116650) message("Time bin=%d new_dti=%d", p->time_bin,
+  // new_dti);
   if (p->gpart != NULL) {
-    p->gpart->ti_begin = p->ti_begin;
-    p->gpart->ti_end = p->ti_end;
+    // p->gpart->ti_begin = p->ti_begin;
+    // p->gpart->ti_end = p->ti_end;
+    p->gpart->time_bin = get_time_bin(new_dti);
   }
 
   /* Get the acceleration */

src/runner.c

@@ -113,6 +113,9 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
 #include "runner_doiact_fft.h"
 #include "runner_doiact_grav.h"
 
+#undef ICHECK
+#define ICHECK -116650
+
 /**
  * @brief Perform source terms
  *
@@ -205,7 +208,7 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
   struct part *restrict parts = c->parts;
   struct xpart *restrict xparts = c->xparts;
   const int count = c->count;
-  const int ti_current = r->e->ti_current;
+  // const int ti_current = r->e->ti_current;
   const struct cooling_function_data *cooling_func = r->e->cooling_func;
   const struct phys_const *constants = r->e->physical_constants;
   const struct UnitSystem *us = r->e->internalUnits;
@@ -227,9 +230,11 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
       struct xpart *restrict xp = &xparts[i];
 
       /* Kick has already updated ti_end, so need to check ti_begin */
-      if (p->ti_begin == ti_current) {
+      // if (p->ti_begin == ti_current) {
+      {
 
-        const double dt = (p->ti_end - p->ti_begin) * timeBase;
+        // const double dt = (p->ti_end - p->ti_begin) * timeBase;
+        const double dt = 1. * timeBase;  // MATTHIEU
 
         cooling_cool_part(constants, us, cooling_func, p, xp, dt);
       }
@@ -906,6 +911,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
 
   const struct engine *e = r->e;
   const double timeBase = e->timeBase;
+  const int ti_current = e->ti_current;
   const int count = c->count;
   const int gcount = c->gcount;
   struct part *restrict parts = c->parts;
@@ -946,15 +952,17 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
           const int new_dti = get_gpart_timestep(gp, e);
 
           /* Now we have a time step, proceed with the kick */
-          kick_gpart(gp, new_dti, timeBase);
+          kick_gpart(gp, new_dti, e->ti_current, timeBase);
 
           /* 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);
+        const integertime_t ti_end =
+            get_integer_time_end(ti_current, gp->time_bin);
+        ti_end_min = min((int)ti_end, ti_end_min);
+        ti_end_max = max((int)ti_end, ti_end_max);
       }
     }
 
@@ -967,9 +975,18 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
       struct part *restrict p = &parts[k];
       struct xpart *restrict xp = &xparts[k];
 
+      /* Minimal time for next end of time-step */
+      integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+      if (p->id == ICHECK)
+        message("Particle in kick ti_end=%lld ti_current=%d", ti_end,
+                e->ti_current);
+
       /* If particle needs to be kicked */
       if (part_is_active(p, e)) {
 
+        if (p->id == ICHECK) message("Particle active in kick");
+
         /* First, finish the force loop */
         hydro_end_force(p);
         if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
@@ -977,17 +994,27 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
         /* Compute the next timestep (hydro condition) */
         const int new_dti = get_part_timestep(p, xp, e);
 
+        if (p->id == ICHECK)
+          message("time_step=%d (%e)", new_dti, new_dti * e->timeBase);
+
         /* Now we have a time step, proceed with the kick */
-        kick_part(p, xp, new_dti, timeBase);
+        kick_part(p, xp, new_dti, e->ti_current, timeBase);
+
+        // if (p->id == ICHECK) printParticle_single(p, xp);
 
         /* Number of updated particles */
         updated++;
         if (p->gpart != NULL) g_updated++;
+
+        ti_end += get_integer_timestep(p->time_bin);
       }
 
-      /* Minimal time for next end of time-step */
-      ti_end_min = min(p->ti_end, ti_end_min);
-      ti_end_max = max(p->ti_end, ti_end_max);
+      if (p->id == ICHECK)
+        message("ti_current = %d dti=%lld ti_end=%lld (%f)", ti_current,
+                get_integer_timestep(p->time_bin), ti_end,
+                ti_end * e->timeBase);
+      ti_end_min = min((int)ti_end, ti_end_min);
+      ti_end_max = max((int)ti_end, ti_end_max);
     }
   }
 

src/statistics.c

@@ -125,7 +125,10 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
     const struct gpart *gp = (p->gpart != NULL) ? gp = p->gpart : NULL;
 
     /* Get useful variables */
-    const float dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
+    const integertime_t ti_begin =
+        get_integer_time_begin(ti_current, p->time_bin);
+    const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+    const float dt = (ti_current - (ti_begin + ti_end) / 2) * timeBase;
     const double x[3] = {p->x[0], p->x[1], p->x[2]};
     float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]};
     if (gp != NULL) {
@@ -206,7 +209,10 @@ void stats_collect_gpart_mapper(void *map_data, int nr_gparts,
     if (gp->id_or_neg_offset < 0) continue;
 
     /* Get useful variables */
-    const float dt = (ti_current - (gp->ti_begin + gp->ti_end) / 2) * timeBase;
+    const integertime_t ti_begin =
+        get_integer_time_begin(ti_current, gp->time_bin);
+    const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
+    const float dt = (ti_current - (ti_begin + ti_end) / 2) * timeBase;
     const double x[3] = {gp->x[0], gp->x[1], gp->x[2]};
     const float v[3] = {gp->v_full[0] + gp->a_grav[0] * dt,
                         gp->v_full[1] + gp->a_grav[1] * dt,

src/timeline.h

@@ -24,8 +24,11 @@
 
 /* Local headers. */
 #include "inline.h"
+#include "intrinsics.h"
 
-typedef long long integertime_t;
+#include <math.h>
+
+typedef unsigned long long integertime_t;
 typedef char timebin_t;
 
 /*! The number of time bins */
@@ -34,6 +37,8 @@ typedef char timebin_t;
 /*! The maximal number of timesteps in a simulation */
 #define max_nr_timesteps (1LL << (num_time_bins + 1))
 
+//#define ceil_div(x ,y) (x + y - 1) / y
+
 /**
  * @brief Returns the integer time interval corresponding to a time bin
  *
@@ -45,6 +50,17 @@ static INLINE integertime_t get_integer_timestep(timebin_t bin) {
   return 1LL << (bin + 1);
 }
 
+/**
+ * @brief Returns the time bin corresponding to a given time_step size.
+ *
+ * Assumes that integertime_t maps to an unsigned long long.
+ */
+static INLINE timebin_t get_time_bin(integertime_t time_step) {
+
+  /* ((int) log_2(time_step)) - 1 */
+  return 62 - intrinsics_clzll(time_step);
+}
+
 /**
  * @brief Returns the physical time interval corresponding to a time bin.
  *
@@ -67,7 +83,10 @@ static INLINE integertime_t get_integer_time_begin(integertime_t ti_current,
                                                    timebin_t bin) {
 
   const integertime_t dti = get_integer_timestep(bin);
-  return ti_current % dti;
+  if (dti == 0)
+    return 0;
+  else
+    return dti * ((ti_current - 1) / dti);
 }
 
 /**
@@ -81,7 +100,10 @@ static INLINE integertime_t get_integer_time_end(integertime_t ti_current,
                                                  timebin_t bin) {
 
   const integertime_t dti = get_integer_timestep(bin);
-  return ti_current % dti + dti;
+  if (dti == 0)
+    return 0;
+  else
+    return dti * ceil((double)ti_current / (double)dti);
 }
 
 /**

src/timestep.h

@@ -23,31 +23,34 @@
 #include "../config.h"
 
 /* Local headers. */
-#include "const.h"
 #include "cooling.h"
 #include "debug.h"
+#include "timeline.h"
+
 /**
  * @brief Compute a valid integer time-step form a given time-step
  *
  * @param new_dt The time-step to convert.
- * @param ti_begin The (integer) start of the previous time-step.
- * @param ti_end The (integer) end of the previous time-step.
+ * @param old_bin The old time bin.
+ * @param ti_current The current time on the integer time-line.
  * @param timeBase_inv The inverse of the system's minimal time-step.
  */
-__attribute__((always_inline)) INLINE static int get_integer_timestep(
-    float new_dt, int ti_begin, int ti_end, double timeBase_inv) {
+__attribute__((always_inline)) INLINE static int make_integer_timestep(
+    float new_dt, timebin_t old_bin, integertime_t ti_current,
+    double timeBase_inv) {
 
   /* Convert to integer time */
-  int new_dti = (int)(new_dt * timeBase_inv);
+  integertime_t new_dti = (integertime_t)(new_dt * timeBase_inv);
 
-  /* Recover the current timestep */
-  const int current_dti = ti_end - ti_begin;
+  /* Current time-step */
+  integertime_t current_dti = get_integer_timestep(old_bin);
+  integertime_t ti_end = get_integer_time_end(ti_current, old_bin);
 
   /* Limit timestep increase */
-  if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
+  if (old_bin > 0) new_dti = min(new_dti, 2 * current_dti);
 
   /* Put this timestep on the time line */
-  int dti_timeline = max_nr_timesteps;
+  integertime_t dti_timeline = max_nr_timesteps;
   while (new_dti < dti_timeline) dti_timeline /= 2;
   new_dti = dti_timeline;
 
@@ -82,8 +85,8 @@ __attribute__((always_inline)) INLINE static int get_gpart_timestep(
   new_dt = max(new_dt, e->dt_min);
 
   /* Convert to integer time */
-  const int new_dti =
-      get_integer_timestep(new_dt, gp->ti_begin, gp->ti_end, e->timeBase_inv);
+  const int new_dti = make_integer_timestep(new_dt, gp->time_bin, e->ti_current,
+                                            e->timeBase_inv);
 
   return new_dti;
 }
@@ -138,8 +141,8 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
   new_dt = max(new_dt, e->dt_min);
 
   /* Convert to integer time */
-  const int new_dti =
-      get_integer_timestep(new_dt, p->ti_begin, p->ti_end, e->timeBase_inv);
+  const int new_dti = make_integer_timestep(new_dt, p->time_bin, e->ti_current,
+                                            e->timeBase_inv);
 
   return new_dti;
 }